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ABSTRACT 


A prime obstacle to the widespread use of adaptive control is 
the degradation of performance and possible instability resulting 
from the presence of unmodeled dynamics. The approach taken is to 
explicitly include the unstructured model uncertainty in the oytput 
error identification algorithm. The order of the compensator is 
successively increased by including newly identified modes. During 
this model building stage, heuristic rules are used to test for con- 
vergence prior to designing new compensators. Additionally, the 
recursive identification algorithm has been extended to multi-input, 
multi-output systems. Enhancements were also made to reduce thp 
computational burden of an algorithm for obtaining minimal state 
space realizations from the inexact, multivariable transfer func- 
tions which result from the identification process. A number of 
potential adaptive control applications for this approach are illus- 
trated using computer simulations. Results indicated that when 
speed of adaptation and plant stability are not critical, the 
proposed schemes converge to enhance system performance. 
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Chapter I 


INTRODUCTION 

A. ADAPTIVE CONTROL BACKGROUND 

Automatic control systems are often required to operate 
physical plants which have a wide range of dynamic properties. 
These changes in dynamic properties could be caused by one of the 
followings 1) alterations in operating environment; 2) structural 
modification of the plant; or 3) failure of one of its compo- 
nents. If the variation of the dynamic properties is sufficiently 
large, a single point design of the control system (no matter how 
robust) may not be able to satisfy the performance specifications. 
Hence, a control system is required that can reconfigure itself to 
provide enhanced performance in the face of these variations. A 
control system is said to "adapt" if its internal structure or 
method of obtaining feedback control is altered in response to these 
changes. 

Examples of problems which would benefit from efficient 
adaptive control systems exist in most areas of engineering. Flex- 
ible flight vehicles may be required to fly over a wide range of 
velocity, height and Mach number which significantly modifies their 
flight characteristics. In the case of jettisonable flight stores, 
rapid changes in dynamic properties occur. Space structures are 
difficult to analyze and test on Earth, have a wide range of chang- 
ing external disturbances and may actually grow significantly during 
deployment and/or construction. Additionally, if for any applica- 
tion some element of the plant or control system fails, a reconfig- 
uration of the control system may be desirable for either safety or 
performance. 


Adaptive control strategies have two subdivisions of effort 
which can be used to distinguish them from nonadaptive systems (See 
fig, (1-1)). The first subdivision is a learning system which 
improves the information about the unknown system variables. 
Another task of this estimation subdivision may be the construction 
of state variables for use by the other subdivision, the control- 
ler. The controller subdivision determines the control inputs to 
the plant based upon the estimated state variables and plant para- 
meters. Adaption occurs when either one or both of the subdivisions 
alter the computation scheme based upon the values of the state 
variables, primary plant parameters or secondary (dependent) plant 
parameters. 

Especially in the case of continuously adapting systems, the 
two subdivisions interrelate. This interrelationship makes conver- 
gence and stability proofs difficult, even in cases where perfect 
modeling exists. The estimation subdivision uses past control 
inputs and past plant output measurements to generate current state 
estimates and parameter estimates. A key feature of this subdivi- 
sion, however, is that the estimation accuracy is heavily dependent 
upon the control inputs. A clever sequence of controls can be used 
to excite specific modes, isolate effects of certain gain parameters 
and regulate signal-to-noise ratios at the sensors. However, con- 
trol inputs that are good for estimation purposes may not be good 
for mission or performance specifications. 

The controller subdivision computes gains based upon the 
current estimate of the state variables and plant parameters. While 
the estimation subdivision utilizes previous control inputs and out- 
put measurements, the controller uses the current values and esti- 
mates of future values to compute current control inputs. If the 
estimation subdivision fails, the controller subdivision will also 
fail. For this reason, most adaptive control research emphasizes 
the importance of accurate and efficient algorithms in the 
estimation subdivision. 
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FUNDAMENTAL SUBDIVISIONS OF ADAPTIVE CONTROL 


U 


Figure 1-1 



Block diagram showing fundamental subdivisions of 
adaptive control. 
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An adaptive system may be defined as one which measures its 
performance relative to some index and modifies its internal para- 
meters to approach a set of optimum values [1]. It has also been 
suggested that an adaptive system is one that is designed from an 
adaptive viewpoint [2]. Clearly it becomes a logical impossibility 
to determine by observation of performance characteristics whether 
or not a control system is adaptive [3]. The very use of feedback 
in control systems reduces the system sensitivities to external 
disturbances and changes in the plant operating characteristics. 
The behavior of a system with feedback tends to be invariant to 
internal changes in itself or in the environment and may legiti- 
mately be called adaptive in the customary usage of that word. 

In the context of control theory, an adaptive system usually 
excludes control system designs in which the state variables are 
measured or estimated and the parameters are assumed to be known. 
Although the distinction between state variables and parameters is 
usually based upon historical perspective, it really is a conveni- 
ence of the control system design engineer and not necessarily an 
objective observation [3]. As an example, consider the following 
process 


x = -px + u (1.1) 

It would be customary to call x the state variable and p the 
parameter. Suppose that p is not a constant and can really be 
modeled as 


p = xv (1.2) 

with v as another input variable (control or external distur- 
bance). Now the designer can approach the problem in two different 
ways. 

If an adaptive system is desired, a controller is designed 
using equation (1.1) only. The parameter p is tracked by an 
identification scheme or estimated through explicit measurements to 
maintain acceptable performance as it varies throughout its range. 
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On the other hand, the designer could implement a nonadaptive system 
considering both (1.1) and (1.2) and explicitly using x and p as 
state variables. It would be easy to evaluate the performance of 
either system, but it would be impossible to determine which 
viewpoint was used by the designer. 

A universal adaptive controller is a fictitious device which is 
the ultimate goal of control designers. It consists of a box with 
inputs and outputs which could be connected to any process. It has 
no prior knowledge of the system it tries to control, but after a 
period of time, it generates an internal model and begins to control 
the plant. It continues to be sensitive to any plant changes or 
variations in the external disturbance field while maintaining good 
performance characteristics. Universal adaptive controllers for 
general plants are assuredly a long way in the future. However, the 
research reported herein is a step in this direction. It represents 
an attempt to identify and build an internal model structure with 
limited prior knowledge for use in adaptive control of a limited 
class of problems. 

In this report, a system will be considered adaptive if it is 
implemented with the following design principles: 1) continuous 
monitoring of system performance; 2) use of a f igure-of-merit for 
decision making; and, 3) adjustment of internal parameters or struc- 
ture to improve the performance. In addition, a distinction will be 
made between parameter adaptive control and adaptive control. Para- 
meter adaptive control is where the control structure or internal 
model order is held constant. In contrast, adaptive control is the 
wider class of problems which requires a change in the order (number 
of state variables) of the internal representation. Most of the 
research in this area has been with parameter adaptive systems. The 
next section outlines some approaches to parameter adaptive control 
and discusses some of the inherent problems when unmodeled dynamics 
are neglected. 


B. CURRENT STATUS OF PARAMETER ADAPTIVE CONTROL 


B-l CURRENT ALGORITHMS 

Adaptive control has received attention from theoreticians and 
practitioners for the past 25 years. Nearly a dozen books and 
hundreds of papers have been devoted to the subject. Most of this 
research has been on parameter adaptive control which maintains the 
order of the modeled system as constant, thereby neglecting 
unmodeled dynamics. Examples of early approaches which rely princi- 
pally upon analog circuit technology can be found in references 
[4 >5] • With the advent of computer technology, new methods of para- 
meter adaptive control were developed, including the following: 
model reference adaptive control, self-tuning regulators, dual- 
control methods, multiple— model adaptive control, and adaptive 
observation. These methods are reviewed below as a framework for 
understanding the problems that are being addressed by this 
research. 

ADAPTIVE OBSERVERS — A logical approach toward adaptive control 
is to use an adaptive observer in conjunction with a constant gain 
matrix multiplying. the estimated states as depicted in 
figure (1-2). Advances have been made which guarantee stability of 
this approach under certain restrictions [6] and [7]. This method 
has been further extended in [8] to include adaptive gain selection 
for the control law and is illustrated in figure (1-3). 

MODEL REFERENCE ADAPTIVE CONTROL— Explicit and implicit model 
reference adaptive control (MRAC) algorithms have proven to be 
extremely useful for controlling plants with wide variations in 
plant parameters [9]. Explicit MRAC is where a reference model is 
specified and an adaptive control algorithm is used to make the 
plant output asymptotically approach the reference model output. 
Implicit MRAC includes the case where the reference model is 
adjusted during the adaptation process as an intermediate step. 
Again the purpose is to drive the output of the plant asymptotically 
toward the model response. Both approaches are illustrated in 
figures (1-4) and (1-5). 
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ADAPTIVE OBSERVER METHOD OF ADAPTIVE CONTROL 



Figure 1-2 — Block diagram of adaptive observer method of 
parameter adaptive control. 


ADAPTIVE OBSERVER WITH ASYMPTOTIC GAIN SELECTION 
FOR ADAPTIVE CONTROL 



Figure 1-3 — Block diagram showing adaptive control via adaptive 
observation with asymptotic controller gain 
selection. 
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EXPLICIT MODEL REFERENCE ADAPTIVE CONTROL 



Figure 1-4 — Fundamental block diagram for explicit model 
reference adaptive control (MRAC) algorithms. 


IMPLICIT MODEL REFERENCE ADAPTIVE CONTROL 



Figure 1-5 — Fundamental block diagram for implicit model 
reference adaptive control (MRAC) algorithms. 
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SELF-TUNING REGULATOR — A structure representing a self-tuning 
regulator (STURE) is shown as figure (1-6). Although there are a 
number of ways to implement this approach [10,11], it usually con- 
sists of the use of output error identification techniques to update 
model parameters. A linear dependence of the control gains upon 
these parameters is typically developed so that the control strategy 
is adaptively updated on-line with the parameters. 

MULTIPLE MODEL ADAPTIVE CONTROL—Multiple model adaptive 
control (MMAC) was recently investigated in references [12,13]. 
Figure (1-7) depicts the basic approach. A number of models are 
assumed and used in the state estimation subdivision. The models 
are evaluated based upon their performance and either 1) the best 
model is used in developing the control law, or 2) a weighted sum of 
the models is used in the control law computation. The first method 
affords more flexibility than pure gain scheduling. The second 
method allows one to hypothetically span the range of model 
parameters with good performance. 

DUAL CONTROL — Dual control methods [14-16] use an iterative 
algorithm which explicitly includes a cost for identification in the 
total cost function. As previously mentioned, identification accur- 
acy is highly dependent upon the control inputs, the state variables 
and the plant outputs. By explicitly including a cost for identifi- 
cation, control energy is expended in order to improve the estimates 
of the states and the plant parameters. Unfortunately, dual control 
schemes are currently computationally burdensome, which has so far 
prevented their use in aerospace applications. 

There are approaches to parameter adaptive control other than 
those illustrated above. However, the information flow is similar 
in all approaches and can be represented by the fundamental scheme 
of adaptive control (fig. ( I— 1 ) ) • In the next subsection the impact 
of unmodeled dynamics upon these adaptive algorithms will be 
considered. 
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SELF TUNING REGULATOR 



Figure 1-6 — Block diagram for self-tuning regulator (STURE). 


MULTIPLE-MODEL ADAPTIVE CONTROL 



Figure 1-7 — Block diagram showing overall structure for 

multiple-model adaptive control (MMAC) algorithm. 
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B-2 CONVERGENCE PROBLEMS OF ADAPTIVE ALGORITHMS 


Although adaptive control has been actively researched for the 
past 25 years, there have been very few cases where it has been used 
in actual control systems • Direct utilization of any particular 

algorithm without fine tuning to the problem is not currently pos- 
sible* Recent research has been directed at trying to understand 
the fundamental problems of parameter adaptive control algorithms 
[17]. For example, in the case of MRAC algorithms problems tend to 
occur in cases where prior knowledge of the plant or operating 

environment is poor [18]. Specifically, MRAC algorithms suffer from 
the following problems: a) generation of high frequency control 

inputs; b) high susceptibility to instability in the presence of 
unmodeled dynamics; and c) poor performance in the presence of 
observation noise. 

Adaptive control algorithms have a number of universal problems 
which need to be resolved. If large reference inputs are used, 

there is a tendency for the adaptation process to react too fast, 

sometimes leading to instability. An implicit adaptive algorithm 

which relies upon internal system identification as part of the 
adaptation process can be proven to converge with global stability 
only when unraodeled dynamics do not exist. In fact, parameter adap- 
tive algorithms characteristically result in high gain/high band- 
width systems where the gain and bandwidth tend to grow in the 
presence of disturbances. Clearly this can be disastrous in terras 
of control spillover for truncated models. 

These unfavorable aspects are a direct consequence of the 
information flow between the estimation subdivision and the control 
subdivision of the adaptation process. Errors in the predicted out- 
put and the commanded output are utilized by both subdivisions to 
simultaneously improve the performance, oftentimes resulting in 
growing bandwidth and gain. Classical solutions such as low pass 
filtering confuse the issue since stability and convergence proofs 
for the currently available adaptive algorithms require known model 
order and a fixed relative order. 
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The analytical research performed in references [17-20] relied 
upon a linearization procedure called M final approach analysis" to 
study adaptive algorithms near the solution. The closed-loop, non- 
linear, time varying equations were linearized for single-input, 
single-output (SISO) first order systems when the system outputs and 
reference model outputs were close. This allowed a study of many of 
the adaptive algorithms for making relative performance compari- 
sons. In addition, an exhaustive set of simulation studies were 
performed on many of the known parameter adaptive algorithms [19]. 
Results of these studies showed that most of the algorithms con- 
verged and were stable under ideal conditions. Some were even 
robust with adverse white noise in the process. However, all of the 
adaptive algorithms studied diverged under the following conditions: 

1. Small bias in controls or sensors 

2. Unmodeled, noncontrollable mode as sensor noise 

3. Unmodeled dynamics in plant 

These conditions are almost always encountered in real system 
implementations. Therefore, current parameter adaptive control 
techniques cannot be used without careful tuning. 

Clearly, the direction of adaptive control research must 
change. If we assume that a controller exists which can handle 
unmodeled dynamics, sensor noise and biases, what should its form 
be? Methods for the estimating parameters of the structured model 
uncertainty are fairly mature; however, emphasis is needed to 
develop parameter identification algorithms for working on what is 
called the unstructured model uncertainty. This implies model 
building (or modification) on-line. A successful adaptive control- 
ler will need to modify or add on to its internal model representa- 
tion whenever the unstructured model uncertainty significantly 
degrades the control system performance. 

It has been argued that a requisite feature of a truly adaptive 
control algorithm which can minimize the risk of failure due to 
unstructured uncertainty is that it must have some level of machine 
intelligence. It should take advantage of computer technology to 
monitor its performance and adjust its characteristics in response 
to unmodeled (but predictable) dynamics or system errors. So a set 
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of heuristic rules may be necessary to insure widespread application 
to a variety of problems with a minimum of specific problem fine 
tuning. The decisions required of an intelligent controller may be 
similar to those an experienced engineer would need to make to opti- 
mize an algorithm to a particular system. 

C. PROPOSED APPROACH FOR ADAPTIVE CONTROL 

The research reported herein is an engineering approach to try 
and solve some of the characteristic problems of adaptive control 
that were mentioned in the previous sections. A fundamental dif- 
ference of this research when compared to previous efforts in adap- 
tive control is that the internal model order of the compensator 
will not necessarily be considered fixed. A systematic approach for 
building a model during real-time processing will be developed. At 
the crux of this development is the extension of an output error 
identification algorithm to explicitly include the unstructured 
model uncertainty. It will be shown that under a number of restric- 
tions, it is possible to implement an algorithm with the capability 
of performing on-line equivalent system identification for a number 
of applications of adaptive control. 

The algorithms are developed assuming that digital implementa- 
tion [21 J for current microprocessors is desired. The implicit 
parameter estimation algorithm that was chosen comes from digital 
signal processing and is referred to as the LMS (Least Mean Square) 
algorithm. The LMS algorithm has a very low computational burden at 
each time step, making real-time, recursive processing in a micro- 
processor possible. The LMS algorithm estimates parameters of a 
z-domain SISO transfer function. In chapter II this is extended to 
multi-input, multi-output (MIMO) systems with an autoregressive 
moving average (ARMA) model format. Since linear quadratic gaussian 
(LQG) multivariable design techniques [22] are used for the esti- 
mator and controller design, it is necessary to transform the MIMO 
z-domain transfer function to a minimal state space form. A compu- 
tationally efficient algorithm for performing the required state 
space realization is developed in chapter IV. 
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Figure (1-8) shows the structure of the model identification 
algorithm. It is similar to the structure proposed in references 
[23,24]. A suboptimal Kalman-Bucy filter, which is suboptimal in 
the sense that prior knowledge of the noise statistics and the 
actual plant are not available, is designed and used to estimate the 
states for possible use by the controller. Simultaneously and in 
parallel, the adaptive algorithm is processing the input and output 
data to identify new parameters or update the ones currently used in 
the model for designing the Kalman-Bucy filter. Periodically the 
asymptotic Kalman-Bucy filter is redesigned either with the updated 
parameters or by adding another mode to the Kalman-Bucy filter model 
(increasing its order). This process is continued until some 
evaluation of system performance is satisfied. 

The closed-loop flow diagram is depicted in figure (1-9). The 
parallel structure between the adaptation process and the controller 
is critical. The separation of the two functions is important as it 
prevents many of the failures mentioned in the previous sections. 
Large inputs, colored noise (possibly due to unmodeled dynamics) and 
biases, which may cause parameter estimation inaccuracies during 
adaptation do not feed directly through to the controller design. 
The properties of LQG optimal compensators in the presence of 
unmodeled disturbances have been extensively studied [25,27]; and 
while system performance may be compromised, the divergence charac- 
teristic of parameter adaptive control algorithms will only occur in 
extreme cases. 

Another advantage of the separation of the implicit model 
identification task from the system control task is that internal 
monitoring is possible. The identification accuracy is principally 
a function of the control inputs used to excite the system. Since a 
computer program is used to monitor system performance and to peri- 
odically update system parameters and/or model order, it can also be 
used to make decisions about how to improve the equivalent system 
identification accuracy, either through adjusting some free para- 
meters in the adaptive identification scheme or through specifying a 
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ADAPTIVE FILTER IN PARALLEL WITH KALMAN FILTER 



Figure 1-8 — Block diagram showing use of adaptive recursive 
identifier for model building in parallel to 
Kalman- Bucy filter. 
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OVERALL BLOCK DIAGRAM 
FOR PROPOSED ADAPTIVE CONTROL SCHEME 



Figure 1-9 — Block diagram for proposed on-line equivalent system 
identification scheme for adaptive control. 




different set of control Inputs. This decision making process is 
equivalent to the first stages of programming heuristic problem 
solving logic. 

A significant contribution of this research is the modification 
to the output error identification algorithm, LMS, to include addi- 
tional states and is called the Incremental Mode LMS (IMLMS) algo- 
rithm. The IMLMS algorithm is developed in chapter II. A number of 
simulations illustrating the parameter identification capabilities 
of the LMS and IMLMS algorithm are presented in chapter III. Com- 
plete examples of equivalent system identification with model 
building and adaptive control are presented in chapter VI. 
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Chapter II 


OUTPUT ERROR IDENTIFICATION - ANALYSIS 
A. INTRODUCTION 

The performance of adaptive control algorithms depends 
primarily on the manner in which the estimation subdivision reduces 
the model uncertainty. Implicit adaptive control techniques use 
parameter identification to formally find an internal model repre- 
sentation of the plant prior to applying a control law. Hence, 
asymptotic stability of the overall algorithm requires that the 
model parameters approach those of the best equivalent system repre- 
sentation of the plant. This is easily understood in the case of 
plants which can be accurately represented by low order, linear, 
time— invariant models but becomes complicated when the plant has 
significant nonlinearities — higher than modeled order or time varia- 
tions. It is unfortunate that these cases of high model uncertainty 
are the very ones which drive the engineer toward adopting adaptive 
control strategies as a means of obtaining or maintaining system 
performance. Adaptive control cannot be an attractive alternative 
to robust control [28] until it can treat both modeled and unmodeled 
plant uncertainties. 

In this research, an implicit adaptive control scheme is 
utilized which requires the on-line identification of plant para- 
meters. Two anticipated applications are aircraft flight control 
(particularly in regions of significant nonlinearities or changing 
parameters) and spacecraft modal damping with pointing control. In 
the case of spacecraft problems, a sample rate in excess of 200 Hz 
has been envisioned [29]. The high sample rate required, coupled 
with the relatively long periods of the rigid body modes, make 
purely batch processing techniques (e.g., least squares or maximum 
likelihood) impractical both from computational burden and data 
storage aspects. This is true especially for microprocessor 
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implementations. For these reasons, a recursive algorithm was 
choosen for the discrete output error identification operation of 
the adaptive control scheme. 

Recursive identification has the advantage of always having an 
updated parameter estimate for use in the controller subdivision. 
However, recursive algorithms are less precise due to the simplifi- 
cations needed to run at high sample rates. Hence, batch processing 
tends to yield more accurate results with fewer samples of data. 
The engineering trade-off is whether or not by the time a batch 
identification scheme has processed the data, a recursive scheme 
could have already converged to an acceptable answer. As previously 
mentioned, the storage and speed requirements of certain applica- 
tions make microprocessor implementation of batch identification 
methods impractical. 

Specifically, the algorithm utilized for output error identifi- 
cation in this research is the LMS (Least Mean Square) Adaptive 
Predictor Filter. It was first introduced by Widrow and Hoff [30] 
and has been further analyzed in references [31-33]. The algorithm 
is very similar to a class of adaptive filters known as SHARF 
(simple hyperstable adaptive recursive filter) [34-36] or sometimes 
HARF (hyperstable adaptive recursive filter) [37-39]. The LMS, 
SHARF and HARF filters are similar in formulation, but differ only 
in the choice of scale factors. References [31-33] analyze the LMS 
algorithm showing a derivation, a study of selecting the magnitude 
of the step size factor, an estimate of the speed adaption, and 
convergence proofs in stationary stochastic environments. Refer- 
ences [37-39] are particularly useful as they also discuss parameter 
convergence in the presence of unmodeled dynamics. 

In this chapter, the output error identification algorithms for 
the adaptive control scheme will be developed. In section II-B the 
single-input, single-output (SISO) form of the LMS algorithm will be 
developed. Section II-C extends the LMS algorithm to multivariable 
(multi-input, multi-output or MIMO) plants. Section II-D studies 
the analytical convergence properties of the output error identifi- 
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cation algorithm. A hybrid batch and recursive version of the 
algorithm is developed in section II-E. A prime contribution of 
this research, an explicit inclusion of the unstructured model 
uncertainty in the output error identification, is included as 
section II-F. 

B. DERIVATION OF SISO OUTPUT ERROR IDENTIFICATION 
In this section the LMS adaptive filter algorithm is derived 
following the development of references [31,32]. It should be noted 
that the notation used in this derivation prevents the extension to 
a multivariable representation, a necessity for the algorithm's use 
in practical stochastic control applications. 

Using the autoregressive moving average (ARMA) representation 
[21], a dynamic system is represented in terms of past measurements, 
y(k), and past control inputs, u(k), as: 


y(k) = x T ( k ) «00 


( 2 . 1 ) 


where 


x T (k) = [y(k-l) ,y(k-2) ... y(k-n) ,u(k-l) ,u(k-2) 

... u(k-n)] (2.2) 

a(k) is the vector of weights multiplying past state measurements 
and controls for obtaining the predicted output, y(k), of output 
y(k) , and n is the order of the model. The error is defined by 
comparing the measured output and the predicted output 


e(k) = y(k) - y(k) (2.3) 

Substituting (2.1) into (2.3), results in 
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(2.4) 


e(k) = y(k) - x T (k)a(k) 


It is desired to minimize the square of the error for all time 
(-oo<k<°°) by finding the optimum set of weights, a*. A logical 
approach would be to update a trial set of weights using a simple 
steepest descent algorithm. 

£Ll+i = SL± ~ ^ » (2.5) 


where the gradient is given by 


V 


9 «i 


l 

j— e 


e 2 (j) 


( 2 . 6 ) 


Since a recursive algorithm is desired, an update for each 
discrete time step is needed. The gradient of the square of the 
error is formulated as a time dependent variable 




k 

l 


e 2 (j) 


(2.7) 


This would still require a great deal of storage and computation for 
each at each step k. The convention is to update the 

at each time step, so adopting the notation 

= a(k) , (2.8) 


and making the following crucial approximation 

V(k) - V(k) = e 2 (k) , (2.9) 

an estimate for the gradient is obtained. That is, instead of 
summing all the past errors using the current value of a(k), just 
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use the current error. The definition of the gradient from (2.6) 
yields 


V(k) = 2e(k) 


8e(k) 

8a(k) 


( 2 . 10 ) 


Using (2. A) 


3e(k) 

3a(k) 


-x(k) 


( 2 . 11 ) 


so (2.10) becomes 


V(k) = -2e(k)x(k) (2.12) 

Substituting (2.12) into (2.5) using (2.6), the LMS adaptive 
filter for updating the weights is obtained 

a(k+l) = a(k) + 2ye(k)x(k) , (2.13) 

or it can be rewritten as 

a(k+l) = a(k) + yx( k )ty( k )“iS. T ( k ) a ( k )] • (2.14) 

This differs from recursive least squares [21] only in the term u, 
which replaces the estimate-error covariance matrix P(k+1). There 
P is updated at each step by 

P - 1 (k+l) = P - 1 (k) + x(k)x T (k), (2.15) 


or 


P ( k+l ) = P(k) - P(k)x(k)[l+x T (k)x(k)]x T (k)P(k) (2.16) 


23 



In effect , the LMS algorithm replaces the step-varying P matrix by 
a constant diagonal matrix pi, where I is the unity matrix. Hence, 
a reasonable estimate for p is 

H“^tr(P), (2.17) 

i.e., the average error variance of the parameters a^. The 
difference between the LMS algorithms and the SHARF algorithms also 
lies in the selection of weighting factors. 

The parameter update from (2.13) can be made very quickly as it 
requires only about 5n operations. Despite its simplicity, it is 
data adaptive and has convergence properties that approach those of 
more cumbersome conventional methods. Some simple convergence 
arguments follow. Take the expected value of (2.12) 

A 

E[V(k)] = E[-2e(k)x-(k)] (2.18) 
Noting that e(k) is a scalar and using (2.4) 


E[V(k)] - -2E[y(k)x^ (k) - jc( k)^ (k)a( k) ] (2.19) 

Defining the stochastic covariances as 

P(k) = y(k)x T (k) 

( 2 . 20 ) 

R(k) = x(k)x T (k) 

It can be shown that the gradient is equal to (e.g., [41]) 

V = Ra - P (2.21) 

Hence 


E[V(k)] = V £ 


( 2 . 22 ) 
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Solving (2.21) for the optimum weight vector gives the well known 
Wiener solution 

a* = R _1 p (2.23) 

A 

Since the mean value of the gradient estimate, V(k) is equal 
to the gradient, V, the estimate must be unbiased* 

The above analysis assumes that the weight vector was held 
constant. References [31,33,41] extend the analysis to include 
variation of weights and show that the Weiner solution is again 
obtained for the optimum weights if the inputs are uncorrelated over 
time and are stochastically stationary. References [31,33] remove 
the stationarity requirement in a formal manner. A discussion of 
the weight vector convergence is given in section II-D. 

This development of the LMS algorithm can be extended to find 
bounds for the estimate of step size factor, p, to be chosen by the 
engineer. The analyses of references [31,32,41] shows that 

0 < * < <2 - 24) 

as a requirement for convergence. Simulation studies in 

reference [32] indicate that the normalized error of misadjustment 
of weight parameters is linearly proportional to the number of 
weights. In addition, a discussion of the speed of adaptation shows 
it to be an exponential function of p, the eigenvalues of R and 
the number of weight parameters. References [31,33] indicate that 
based upon experience, a “good" value for p for design purposes is 
around .001 plus average input signal power. 

Reference [33] is significant in that it discusses the LMS 
algorithm in terras of a wider class of problems and allows that 
adequate performance is possible under dependent random environ- 
ments, slowly changing parameters and in the presence of unmodeled 
dynamics. A prime contribution of the literature for SHARF and HARF 
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algorithms is the explicit inclusion of the recursive identification 
scheme in an adaptive control function. Analyses in these refer- 
ences [34-39] are concerned with proving global asymptotic conver- 
gence of the parameters, studying the impact of unmodeled dynamics 
and quantifying the sufficient excitation requirements for parameter 
convergence. Some of these issues will be discussed in section D, 
The next section develops a MIMO formulation of the LMS adaptive 
filter. 


C. MIMO FORMULATION OF IMS ALGORITHM 

Most practical applications of adaptive control involve MIMO 
systems, that is, systems with more than one control input and more 
than one output. In this section the LMS algorithm is extended to 
MIMO systems. The vector notation utilized in the previous section 
will not accomodate multi-output formulations, so a tensor notation 
is used to maintain the same AR (autoregressive) coefficients for 
each output. In addition, a distinction shall be made between the 
terms multiplying the measurements and the controls. 

Let 


y^(k) = £th measurement at time k 


e^Ck) = jth input at time k 

Assume that the £th measurement can be predicted as 
/\ 

y^(k) = a^Ck-l) + a 2 y 2^ k ~ 2 ^ + ... + a^Ck-n) 

1) ^2 j £^2 2 ) + ••• + P n j£ e j (k-n) (2.25) 

the notation can be condensed by making the following definitions 
consistent with an ARMA model: 
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x 1£ (k) = y £ (k-i) 


(2.26) 


^.(k) = ej (k-i) (2.27) 


Tensor 

notation with the 

summation 

summation limits will be 

used: 


n 

Y 

m 

Y 


L 

i=l 

l 

J-l 

where 

n is the model 

order, m 

and r 

is the number of 

outputs* 


y £ (k) = a i (k)x i£ (k) 


convention implying the following 
r 

l 

£=1 

is the number of control inputs 
Now (2.25) can be rewritten 

• e ij£ (k)u i;j (k) . (2.28) 


The error of prediction becomes 


e £ (k) = y £ (k) - y £ (k) (2.29) 

or 

e £ (k) = y £ (k) - x A (k)a i£ (k) - 8^ £ (k)u 1 j (k) . (2.30) 

As in the previous section it is desired to find an approximate 
gradient of the square of the error with respect to the weights, 
ctj/k) and £,(k) . 

So assuming that the gradient can be approximated by using the 
current weights and errors, 

(e £ e P = 2 e J (2.31) 

= 2 e £ 3^ ( y “ a i x i£ " B ij£ u ij )e £ (2 ‘ 32) 


including time dependence 
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a ct ± C k > ( e ^ k ) e ^ k )) 2e £ (k)x i£ (k) , 


(2.33) 


and for the control influence terms 




( G o e ») 2e. Can ) ( 




ij* 


u i E t? ' 2E t (y ‘ Vu ‘ VV 


( G ff T (k)e (k)) = -2e (k)u..(k) . 


8 B ij£ 1 * 


£' 7 ij 


(2.34) 


(2.35) 


(2.36) 


So again using a steepest descent approach, equation (2.5) can be 
used to give 


3(e T (k)e (k)) 
“i a+1) ' a l (k) - U 3^0) 

a i (k+l) = a i (k) + 2pe^(k)x iJ ^(k) 

and 


(2.37) 


(2.38) 


P ij£ (k+1) = P ij£ (k) " v 


9 ^ E £ T(k)e £ (k) ) 


(2.39) 


P ij£ (k+1) = P ij£ (k) + 2pe Jl (k)u 1 j(k) . (2.40) 

The parameter update equations which comprise the MIMO LMS 
algorithm are (2.38) and (2.40). Previous formulations implied a 
separate set of c^'s for each output, where as now there is a 
single set of parameters constrained to remain valid for all outputs 
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while maintaining the desirable convergence properties of the SISO 
LMS algorithm. This permits a simpler, multivariable ARMA formula- 
tion for the plant. The a and 8 parameters are merely the deno- 
minator and numerator terms of a discrete, multivariable transfer 
function. Chapter IV develops an efficient way to compute an 
approximate minimal state space realization from the experimentally 
determined transfer functions. The next section discusses conver- 
gence of the LMS algorithm, applicable to both SISO and MIMO 
systems. 


D« NOTES ON OUTPUT ERROR IDENTIFICATION CONVERGENCE 

In this section, the notation of the previous section will be 
used to develop a first order example to illustrate many of the 
convergence properties of the LMS algorithm. The plant is assumed 
to be known and the convergence of the identifier algorithm is to be 
investigated. 


Plant: 


x(k+l) = ax(k) + bu(k) 


(2.41) 


Estimator: 


A 

x(k+l) = a(k)x(k) + 8(k)u(k) 


(2.42) 


Error: c(k) = x(k) - x(k) (2.43) 

e(k) = (a-a(i-l ) )x( i-1 )+(b-8(i-l ) )u(k-l ) (2.44) 

LMS Identifier: a(k+l) = a(k) + 2 y{ [ a-a(k) ] x(k) 

+ [b-p(k)]u(k) }x(k+l) } (2.45) 

B(k+1) = 3(k) + 2p{ [a-a(k)]x(k) 

+ [b-8(k)]u(k) }x(k+l) } (2.46) 

Equations (2.42) and (2.43) become 


a(k+l ) 
3(k+l ) 


a(k) 

3(k) 


+ 2p 


x(k)x(k+l ) 
x(k)u(k+l ) 


u(k)x(k+l ) 
u(k)u(k+l) 


a-a(k) 

b-B(k) 

- 


(2.47) 
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From (2.47) it is easy to see the requirement for sufficient 
excitation. [c.f . ,36 ,39 ,40 ,42] . If steady state conditions are 
reached such that x(k) = x(k+l) or u(k) = u(k+l) no improvement 
is made in the parameters no matter how large the parameter errors, 
a-a(k) and b-3(k), are. Control inputs are required to improve 6. 

If a step input of magnitude u Q is applied at k=0 with 
x(0)=0, it follows that 


. bu 

x(k) = (1-a ) -7— 2- 
1— a 


(2.48) 


indicating that any stable parameterization of a and 8 with a 
zero frequency gain of b/(l~a) will show x^-x [42]. Clearly, a 
and 8 may not have converged to a and b. 

Combining (2.47) and (2.48) near k=0 , we get 


a(k+l) 


'a(k)" 


0 0 

a - a(k) 


= 


+ 2u 

o 


3(k+l) 


3(k) 


0 u 2 

o J 

b - 3(k) 


(2.49) 


whereas for any k, 


“ 



r 

a(k+l ) 

= 

a(k) 

+ 2v 

3 ( k+1 ) 


B(k) 



kv2 


bu 
o 

1-a 

2 


, kv bu o 

(1 ' a > PT 


, abu 
/, k^ o 

°- a > T7 


2 n 


u 


a-a(k) 

b-0(k) 


(2.50) 


and for k+°° 


— — 


- 

r 

a(k+l) 

= 

a(k) 

+ 2 p 

3 (k+1 ) 


3(k) 



bu 


bu 
o 

1-a 

2 


1-a 


abu 


2n 


1-a 

2 

u 


a-a(k) 

b-0(k) 


(2.51) 


For small k the only improvement is in 0 and it can have a 
relatively large magnitude depending upon the sizes of u and Uq. 
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As x(k) moves away from 0, corrections in a begin. It is 
possible, however, for a to initially move in the wrong direction 
depending upon u Q , b-3(k) and p. 

Considering equations (2.49) and (2.50), an input strategy 
suggests itself for obtaining sufficient excitation for parameter 
convergence. A step input provides immediate improvement in p. 
Experience, such as the results in chapter III, indicates that the 
first few samples after the step impulse provide the majority of the 
improvement. As the system begins to respond to a step input, the 
improvement in a grows in magnitude; but the improvement decreases 
at an exponential rate, albeit slower than the improvement in 3. 
If only a single step pulse is applied through the control, the 
parameters will most likely not converge to an accurate answer. 

Multiple step pulses are required for sufficient excitation. 
This is consistent with the general identification research for 
optimal inputs of references [43-45] and the specific LMS results of 
[33,36,37,39]. The pulses are needed to improve the control influ- 
ence terms, 3, but improvement in the denominator terms, a, comes 
only after waiting for the system to respond. So the designer is 
confronted with a conflict when choosing input signals to enhance 
the parameter estimate. Frequent pulses would be advantageous for 
3, but may not give a enough time to begin its improvement. 

Conversely, allowing the system to respond to infrequent pulses 
would yield an ideal convergence environment for a but would be 

insufficient for 3. 

In chapter IV it will be shown that equation (2.42) has a 
simple representation as a discrete (z-domain) transfer function and 
that a corresponds to the denominator term and 3 to the numera- 
tor term. The need for several rapidly changing control inputs to 
estimate the numerator terms is typical of identification problems. 
In contrast, a low frequency externally applied signal (e.g., a 

dither signal) is advantageous for identifying the denominator 

terms. Although white process noise is capable of providing a para- 
meter convergence for a, studies [31,32] indicate that a control- 

induced, signal- to-noise ratio of 2 or more will start to yield the 
exponential convergence rates indicated in equation (2.50). 
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It is apparent that the rate of convergence is dependent upon 
the relative sizes of u and the mean values of x and u. The 
implication is that p should be sized to handle the expected 

values of x and u. An alternative approach, however, may be to 

scale the outputs of the state variable representation so as to 
maintain expected values for x and u near unity. This later 
case is the approach adopted in this research. 

E. BATCH LEAST SQUARES NUMERATOR DERIVATION 

In the previous section it was observed that the improvement in 
the numerator term of the discrete transfer function representation 
of the LMS adaptive algorithm comes primarily within the first few 
samples following a pulse. Therefore, many pulses are required to 
obtain convergence of the numerator terms. However, it is advanta- 
geous to wait on the order of a time constant between pulses so the 
denominator terms can approach their true value. The result is that 
the denominator terms converge more rapidly than the numerator 
terms. 

Since for most, adaptive control applications the speed of 
adaptation is critical, a hybrid batch and recursive identification 
scheme is proposed. Once the denominator terms have been identified 
it becomes a relatively simple matter to use a batch least squares 
algorithm to identify the numerator terms. This operation is tanta- 
mount to finding the zero frequency gain of the system; hence, data 
from the last two pulses are probably all that are necessary for 
sufficient accuracy. It has already been argued that collecting 
data for several system time constants and applying a computation- 
ally burdensome batch identification scheme (e.g., maximum likeli- 
hood techniques [46]) for finding a best model order and/or para- 
meter estimates is inappropriate for microprocessor applications. 
This is true from both a data storage and a time for computation 
standpoint. 

Using a batch scheme for the numerator terms has some potential 
advantages. A prime problem with least squares batch algorithms are 
the inconsistent results obtained with model truncation. This is 
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not as significant a problem for the LMS algorithm, as will be 
illustrated in section III-C. Also, batch data that use only the 
last two pulses are all that is required to be stored for the 
accurate estimation of the control influence terms. Since the 
parameters for the denominator terms from the LMS algorithm are 
available at each sample, the innovation sequence can be recursively 
modified to contain only the numerator content, making the solution 
for the numerator terms a simple computation prior to obtaining a 
state space realization. 

Using an ARMA model representation [21], 

* n n' 

jKk) = - l a £(k-i) + l B. _u(k-i) + e(k) (2.52) 

i=l i=l 1 


Note that is summed up to n 1 , the assumed order of the 

numerator terms, where n f _< n. has the dimensions of r x m, 

where r is the number of measurements and m is the number of 
control inputs. So rewriting (2.52) in tensor notation from 

section II-B 


y £ (k) = -o^Ck-i) + B Uj(fc-i) + e ; (k) (2.53) 

Assume that the are known and that only the part of error 
associated with not knowing the Bij£ is needed. The part of the 
output attributable to the control input u can be written 

y£(k) = y^(k) - a i y Ji (k-i) (2.54) 


Hence, we define a parameter vector, with dimension (nm x 1), 


Q l ^lU P 2U*** B nU e i2A“ ,B nm]l^ T ’ 


(2.55) 


a control vector of dimension (n'm x 1) 
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> 


(2.56) 


mT 


u(k) = tu 1 (k-l)u 1 (k-2)...u 1 (k-n , )u2(k-l)...u m (k-n')] 


and a data vector with N samples with dimension (N-n+lxl) is 
defined as. 


Y^(N) = [y£(n)y^(n+l)...y^(N)J 


(2.57) 


By forming the following matrices. 


U(N) = [u(n), ju(n+l), ...u(N)] 


e £ (N, 6) = [e Jl (n),e Jl (n+l)...e Jl (N)] 


(2.58) 

(2.59) 


The N-n+1 equation errors can be written in matrix notation as 


Y|(N) = U(N)0£ + g a (N,0) 


(2.60) 


The cost function is the square of the errors and is defined by 


W = li (N >V-£ (N ’ 0 P • 


(2.61) 


The normal equations for this case are 


U T WU0£ = u t wy 


(2.62) 


yielding the conventional weighted least squares solution 


0 £ = (u t wu)" 1 u t wy ji . 


(2.63) 
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The weighting matrix, W, is assumed to be the identity matrix 
for the remaining analysis. All that is necessary for this computa- 
tion is to save Nn’m control inputs in the vector U(N) and £N 
samples of the control contribution to output, y f £, in vectors 
y f £(N). It should be noted that y 1 £ is a byproduct of the LMS 
algorithm computation. The a ^ used in the LMS algorithm should 
be satisfactory, provided that enough time constants have elapsed 
and little adaptation is taking place. Finally, it should be 
observed that the 0£ need to be identified for each output, £. 
This comes from the formulation of (2.61) which was adopted for 
saving computer storage space. 

F, INCREMENTAL MODE LMS ALGORITHM DERIVATION 

LMS adaptive filters perform well in the presence of higher 
frequency, unmodeled modes since they act as low-pass filters. This 
capability of the LMS algorithm is studied analytically in refer- 
ences [33,39,40,47] and by simulation in section III-C. It is 
clear, however, that there are times when it would be advantageous 
to add modes to the assumed model structure. Furthermore, if some 
prior knowledge of the plant is available, it would be beneficial to 
include it in the identification process. Chapter I pointed out the 
intrinsic problems of unstructured model uncertainty in terms of 
stability for adaptive control algorithms. In this section an 
approach for adding modes to the model structure will be developed 
by distinguishing between the known (or already identified) part of 
the dynamics and the unknown (or incremental) part. 

One reason that most parameter convergence studies for output 
error identification algorithms are limited to low order is because 
of the inherent numerical inaccuracies that predominate as the order 
is increased. Equation (2.64) is the discrete (z-domain) transfer 
function of equation (2.25) 




+ e 2U z ~ 2 4 — B nn z ~" 


, -1 “2 ^ -n 

1 - a.z - a.z + ... a z 
12 n 


(2.64) 
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The effect of a small error in on the location of the poles 

increases with n. 

A natural approach would be to find the optimum set of a* 
using the LMS algorithm for n, then assume that the order will be 
increased to n+p, and compute a starting set of n+p values of a’ 
using the n values of a* obtained from the LMS algorithm and 
initial guess for p values of a® for the new modes to be 
identified. The indicated product is given by 

l-ct.'z -1 z- (n+p) = (l-a*z -1 +...+a*z“ n ) 

1 n+p I n 

• (l-a?z“ 1 +...+a°z -p ) (2.65) 

v 1 P 

The LMS algorithm can then be used to identify n+p values of a ’ , 
the new coefficients of the larger model. However, this does not 
solve the numerical accuracy problem of higher order systems posed 
above . 

Instead, a method whereby only the new modes are adapted and 
the previously identified modes are held constant is desired. This 
would tend to minimize the impact of numerical inaccuracies of the 
adaptation process upon the estimated system dynamics. It also 
provides a way to include the known dynamics into the identifica- 
tion. The algorithm that is derived below is termed the Incremental 
Mode LMS (IMLMS) algorithm since it provides a way to add incre- 
mental modes to the assumed model. 

We cast the form of the discrete transfer function of (2.64) 
into one which explicitly distributes dynamics into two parts: (1) 

the nonvarying or known part and (2) the incremental or unknown 
part. 


H(z) 


(1 - 


n+p 
l B.z 
i=l 


-i 


l a z i )(l - ^ a z f ) 

i=l f=l 


( 2 . 66 ) 
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where is a rxm matrix for a multivariable system and is 

identified each time, 


e ill B il2 


ilm 


P 21 3 22 e i2m 


P irl e ir2 


irm 


(2.67) 


The a± are assumed constant and only the ctf, the coeffi- 
cients of the incremental modes are identified. The new system 
order is taken as n+p. Multiplying the denominator out in (2.66), 
the following expression is obtained 


H(z) = 


n+p 

r B iZ 

i=l 


-i 


S -i ? -f ? ? -i 

1 - l a z - I a z + l l a ct z 

i=l 1 f=l 1 i=l f=l 


( 2 . 68 ) 


-i-f 


Writing (2.68) in finite difference form [21]: 

y^(k) = a i y £ (k-i) + a f y £ (k-f) - a i a f y £ (k-i-f ) + 8 i£ jUj(k-i) (2.69) 

The implied summation limits of i are 1 to n, of f are 1 to p, 
of % are 1 to r, and of j are 1 to m. Using equations (2.28), 
(2.37) and (2.29) with the above notation a form for the parameter 
update is 


dc, 

a f (k+l) = a f (k) - 2yc £ (k) (k) 


Finding the first partial of the error, e£(k) 



<k) = 


(2.70) 


(2.71) 
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(2.72) 


8e £ 

9aJ ^ = ” a i y £ (k-i-f) 

So the unknown parameter update equation becomes 

<x f (k+l) = a f (k) + 2pe Jl (k)[y Jl (k-f) - a^k-i-f )] (2.73) 

The parameter update equations for the numerator terms, 
remained unchanged as 

0 iij (k+1 > = B i£j^ k ) + Zue^.U) . (2.40) 

The LMS algorithm has been extended to include only the 
adaptation of the unknown modes of arbitrary order. No assumptions 
used in proving asymptotic convergence of the basic LMS algorithm 
have been violated; hence, we expect the IMLMS algorithm to exhibit 
the same asymptotic robustness as the unmodified LMS algorithm. The 
derivation of the IMLMS algorithm was accomplished in the MIMO form- 
ulation of section II- C. Since the numerator update equations are 
unchanged, the same batch least squares identification of the numer- 
ator terms is still possible using the development in section II-E. 

While it is anticipated that the IMLMS algorithm will exhibit 
robust performance under nonideal conditions (sensor noise, proces- 
sing noise or unmodeled dynamics), it will be similar to other iden- 
tification algorithms in that global convergence proofs will not be 
possible in such circumstances. The inclusion of the IMLMS algo- 
rithm in an adaptive control application will need heuristic tests 
to check for divergence prior to utilization of the newly identified 
model. Furthermore, there is a danger of incremental fitting of 
multiple modes which have eigenvalues with nearly equal magnitudes. 
Section III— C will illustrate that spectral separation is necessary 
for good performance in the presence of unmodeled dynamics, enabling 
the IMLMS algorithm to distinguish between modes. 
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The IMLMS algorithm should be especially suitable for aerospace 
applications. Flight vehicles often have a single mode which is 
relatively easy to estimate a priori, but have additional modes 
which are susceptible to wide variations (e.g., phugoid and short 
period modes or roll and Dutch roll modes). Spacecraft typically 
have rigid body modes which are easy to predict, but the interaction 
with the flexible body modes is difficult to estimate or even mea- 
sure with experiments prior to deployment. In both types of vehi- 
cles, the use of an algorithm which holds part of the model constant 
while identifying the incremental part would be advantageous com- 
pared to approaches which identify the entire model for a given 
order. 

In this section, the IMLMS algorithm was presented as an 
effective way to identify explicitly part of a model while holding 
the dynamics of another part fixed. There are many applications 
where this may be a useful approach. In the next chapter, the LMS 
algorithm parameter convergence properties will be examined by simu- 
lation. First, the standard LMS algorithm performance will be 
studied from the viewpoint of a designer. Then, the impact of using 
the multivariable formulation will be ascertained; the influence of 
unmodeled modes will be evaluated; and, the robustness of the IMLMS 
algorithm will be demonstrated. 
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Chapter III 


OUTPUT ERROR IDENTIFICATION - SIMULATION STUDIES 

In this chapter the results of simulation studies investigating 
the use of the LMS algorithm as a parameter estimator are pre- 
sented, These studies are useful in that they can be used to help 
guide the designer in choosing the free parameters of the LMS algo- 
rithm for particular applications. A description of the techniques 
utilized for the computer simulations in this report is given in 
appendix I. The Fortran code for the LMS and IMLMS algorithms is 
given in appendix II. 

A. SISO CHARACTERISTICS OF LMS ALGORITHM 

A-l First Order Demonstration of LMS Algorithm 

In this subsection a first order model will be studied to 
illustrate the parameter estimation convergence rate of the LMS 
algorithm. Several control input strategies are used to excite the 
system. 

The roll mode of an airplane is an example of a first order 
system which is spectrally separated from the remaining lateral- 
directional dynamics (much faster) and may require on-line identifi- 
cation since the damping is a function of angle- of- at tack and 
dynamic pressure. The linear, time-invariant plant for this example 
is modeled by 


ax + bu + w 


(3.1) 


where the baseline value of a is .5 sec"^ and b is 2 
sec”l . Note that the time scale has been transformed for 
convenience as values for a of 10 sec”* are more typical for 

aircraft. The LMS filter equations become 
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x(k+l) = a(k)x(k) + f?(k)u(k) (3.2) 

e(k) = x(k) - x(k) (3.3) 

a(k+l) = a(k) + 2yc(k)x(k) (3.4) 

3(k+l) = 3(k) + 2ye(k)u(k) (3.5) 


A simulation was performed for 50 seconds at a sample rate of 
5 Hz with the inputs scheduled by table 3.1. The step size factor, 
y, was chosen to be .005. A step input was applied at t = 0; from 
10 to 20 seconds a sine wave was superimposed; from 20 to 30 seconds 
a square wave was superimposed; from 30 to 40 seconds the control 
consists of discrete random inputs; and, from 40 to 50 seconds dis- 
Crete process noise is added. Time histories of x and x are 
presented in figure (III-l) and the control input time history is 
shown in figure (III-2). 

The time histories of the estimated parameters are depicted in 
figures (III-3) and (III-4). It is possible to verify some of the 
analytical results of the previous chapter. The denominator coeffi- 
cient, a, moves in the wrong direction initially while the numerator 
coefficient, 3, makes some improvement after a step input. Then a 
begins to move in the correct direction while 3 reaches a steady 
state. Little real benefit comes from the sine wave input while the 
square wave input yields the most improvement in the parameters. 
Random inputs have little impact for this case as signal-to-noise 
ratios of approximately 10 were used. A good rule of thumb is to 
expect parameter convergences with the LMS filter for signal— to— 
process— noise ratios of 2 or more and signal— to— measurement— noise 
ratios of 5 or more [33]. 

At the end of this simulation, the adapted parameters had the 
values shown in table 3.2. Fairly good agreement between the theo- 
retical and experimental coefficients were obtained indicating that 
the identification inputs and other problem parameters were well 
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Table 3«1 — Input schedule for first order demonstra- 
tion of LMS algorithm. 


TIME 

CONTROL INPUT 

NOISE INPUT 

■SB 

u 

w 

0-10 

u 

con 

0 

10 - 20 

u + u sin(fl-lO) 

con amp 

0 

20 - 30 

U + u sgn(sin(fi-10)) 

con amp 

0 

30 - 40 

Tl(0,.l) 

0 

40 - 50 

0 

Tl(0,.l) 


Table 3.2 — Comparison of theoretical and experimentally 
determined parameters for first order 
demonstration of ms algorithm. 


PARAMETERS 

THEORETICAL 

EXPERIMENTAL 

a 

.9048 

.9040 

e 

.3807 

.3857 
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RESPONSE TO CONTROL INPUTS AND PROCESS NOISE 



Figure III— 1 — Time histories of state and state derivative for 
various control and noise inputs. 


CONTROL INPUT TIME HISTORY 



Figure III— 2 — — Time history of control input for first order 
demonstration of LMS filter. 
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ADAPTION OF DENOMINATOR FILTER COEFFICIENT 
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Figure II 1-3 — Time history of LMS filter denominator coefficient 
for first order demonstration. 


ADAPTION OF NUMERATOR FILTER COEFFICIENT 


NUMERATOR 

COEFFICIENT 



Figure III-4 — Time history of LMS filter numerator coefficient 
first order demonstration. 
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matched. In the next subsection, the impact of choosing these free 
parameters will be illustrated. 


A— 2 Study of Important Identification Parameters 

The first order system of the previous subsection is studied 
further in this subsection. The importance of the problem para- 
meters to the identification accuracy is determined by systematic 
variation of each parameter. The parameters considered and their 
nominal values are listed in table 3.3. 

A normalized total parameter error, ep, is defined by 
comparing the estimated parameter value, p^, with the actual 
parameter, p^, and summing for all parameters as follows 


p i “ Pi 

p. = -= =■ x 100% 


(3.6) 


e 

P 



(3.7) 


In this case the parameter vector, j), includes only two items 


p i = 


[o,6] T 


(3.8) 


The system is excited by a square wave for 40 seconds and then 
allowed to settle for 10 seconds. At the end of 50 seconds the 
estimated parameters, a and 8, are compared with their theoretical 
values and the total percent error, ep, is computed. This was 
done for several cases by fixing the values of the problem para- 
meters to their nominal values and varying each one at a time. 

The influence of varying p, the step size factor of the LMS 
algorithm, is depicted in figure (III-5). It shows the expected 
trends. Above a certain value (p “ .8 for this case) the algorithm 
diverges. Two minimum values occur, one for best fit of a and the 
other for 8. This indicates that possibly a different p may be 

desired for each coefficient. It is common to have a different p 
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IMPACT OF STEP SIZE FACTOR UPON CONVERGENCE OF LMS FILTER 



STEP SIZE FACTOR, p 


Figure III-5 — Normalized parameter error versus step size factor 
for first order demonstration of LMS filter. 


IMPACT OF AMPLITUDE OF EXCITING INPUT UPON IDENTIFICATION 



0 1 2 3 4 5 6 


AMPLITUDE OF SQUAREWAVE INPUT. U AMD 

AMP 

Figure III-6 — Normalized parameter errors versus amplitude of 

square wave input for first order demonstration of 
LMS filter. 
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for denominator terms and numerator terms for the SHARF and HARF 

algorithms [36]. If the value of p is too small, the speed of 
adaption is slowed, preventing satisfactory parameter convergence. 
The choice of P is highly problem dependent and is difficult to 
select optimally a priori. 

The amplitude of square wave input has a similar double minimum 
in e p as shown in figure (III-6). If the amplitude of input is 

too large, the algorithm diverges. This is equivalent to too large 
of a stepsize factor, as 0 over-corrects in equation (3.5). If 
the amplitude of the input is too small, sufficient information for 
0 is not available and the system response, x, is too small for the 
nominal value of p. 

The parameter estimation error is very sensitive to the 
frequency of the square wave input as depicted in figure (HI-7). 
The optimum frequency is at .5 rad/sec, which is the system time 
constant. Although a square wave has a wide spectrum of sine wave 
frequencies, it may be advantageous to select square wave inputs 

with varying frequencies. Reference [43] indicates that pulsing the 
controls hard against the stops may be optimum for many linear sys- 
tems, with the switch time for optimum identification a function of 
system parameters. The data in figure (III-7) seem to be consistent 
with this result. 

Although the data have inconsistent variations between 5 and 
30 Hz, figure (III-8) verifies what is expected about system per- 

formance with respect to sample rate. Sampling too slowly results 
in estimation problems due to aliasing effects with respect to the 
input square wave frequency. Sampling too fast runs into numerical 
problems for the nominal values of problem parameters and the accu- 
racy is somewhat reduced. However, it is wiser to sample too fast 
than too slow as the penalties are less and divergence is not a 
problem. 

The plant parameters, a and b, were also varied and the 
estimated parameter errors are shown in figures (III-9) and 
(III-10). Changing a had a much bigger impact than changing b. 

If a is increased in magnitude beyond -1.5, the error becomes 
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IMPACT OF FREQUENCY OF SQUARE WAVE INPUT UPON 
IDENTIFICATION ACCURACY 



Figure III-7 — Normalized parameter error versus frequency of 

square wave input for first order demonstration of 
LMS filter. 


IMPACT OF SAMPLE RATE UPON IDENTIFICATION ACCURACY 
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IMPACT OF SYSTEM DENOMINATOR COEFFICIENT UPON 
IDENTIFICATION ACCURACY 



Figure III-9 — Normalized parameter errors versus plant denomina- 
tor coefficient for first order demonstration of 
LMS filter. 


IMPACT OF PLANT CONTROL POWER COEFFICIENT UPON 
IDENTIFICATION ACCURACY 


NORMALIZED 

PARAMETER 

ERROR. 


% 



Figure III-10 — Normalized parameter errors versus plant control 
influence coefficient for first order demonstra- 
tion of LMS filter. 
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significant, while fairly large changes in b have a much smaller 
impact upon the identification accuracy. 

All of the trends mentioned in this subsection illustrate the 
importance of the input (noise or control) signal on identifiabil- 
ity. The designer must select input signals to sufficiently excite 
the system. In the absence of prior knowledge of system character- 
istics, it may be desirable to use digital implementation of learn- 
ing concepts to vary the input signals. Such heuristic approaches 
are currently under consideration for adaptive control applications. 

The studies of this subsection also show that the step size 
factor, y, is of paramount importance for obtaining good performance 
of the LMS algorithm. It was observed that there are optimum 
values and that if the values are too large, divergence may occur. 
The step size factor may also be a prime candidate for systematic 
variations by computer learning logic. Separate p's for each term 
being adapted may be appropriate. 

A-3 Study of a Second Order System 

In this subsection a second order system will be considered. 
The LMS algorithm implementation will have four terms to be identi- 
fied: two denominator system terms; and, two numerator control 
influence terms. Reference [32] points out that the parameter error 
is a linear function of the number of parameters, so it is expected 
that the accuracy of parameter estimation may be reduced relative to 
the first order studies considered previously. 

A typical second order system that is of concern to the 
airplane control system designer is the short period mode. It is 
spectrally separated from the phugoid mode and its damping ratio 
varies significantly with aircraft geometric characteristics and 
flight condition. Such a system could be generically modeled in 
modal coordinates as 


x = Ax + Bu 


(3.9) 


y = Hx 


(3.10) 
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where 


Kmp Sgn 
U 1o 

If |sin(fit)| 
If |sin(ftt)| 

> .707 
< .707 

(3.11) 

A = 

0 

-a, 2 

1 

— 2Cto 



(3.12) 

B 

= to 

U T 



(3.13) 

H 

= [1 

0] 



(3.14) 


The control input (3.11) was a pulsed wave alternating between 
equal segments of u am p, 0, and _ u am p. The nominal values for 
the design parameters were: y=.005, u amp =.5 sec -2 , fi=l rad/sec, 
o) s =20 Hz and the length of the simulation was for 200 seconds. 
The plant natural frequency, w, was 1 rad/ sec. Only the damping 
factor, C, was varied. Its impact upon average parameter estimation 
accuracy is computed based upon the theoretical estimates for the 
parameters. The results are plotted in figure (III-11). 

The lowest error for the denominator terms of the discrete 
transfer function occurred for a damping ratio of .5. The accuracy 
of the numerator terms increased with the damping ratio and was less 
than the denominator terms when the system was overdamped with £ 
greater than 1.3. These results are consistent with the previous 
observations that the input signal is of prime importance for 
obtaining sufficient excitation for identification. The signal and 
the plant need to be matched for optimum estimation performance. 

The denominator terms converge to accurate results more readily 
than the numerator terms when C is less than 1 . This is a charac- 
teristic of using the LMS algorithm for ARMA type model implementa- 
tions. As previously discussed, the improvements in the estimates 
of the numerator terms, the S's, comes within the first few samples 
of the step input. In contrast, the denominator terms, the a's, 
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IMPACT OF DAMPING RATIO ON LMS IDENTIFICATION 
ACCURACY FOR SECOND ORDER SYSTEM 

12.5 

10.0 

AVERAGE 
PARAMETER 7.5 
ESTIMATION 
ERROR. 

V 5 -° 

% 

2.5 


0 

DAMPING RATIO. £ 


Figure III-ll Average normalized parameter errors versus damping 
ratio for second order system. 
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make the majority of the improvement after the system has had time 
to react. This suggested a hybrid technique using both the recur- 
sive and the batch least squares identification schemes. This 
approach to improve the estimation of the control influence 
parameters was developed in section II-E. 

A-4 Phase/Gain Evaluation of LMS Filter 

Phase and gain plots of the LMS filter were obtained as a 
function of frequency. Such Bode plots are potentially useful for 
evaluating the effect of sampling rate on filter performance. 
These approaches may also be useful for determining the potential 
application of the LMS adaptive algorithm to nonlinear problems. In 
this study, the frequencies of the input signal and the sample rate 
were varied. If describing functions for nonlinear applications 
are desired, the magnitude of the input signal should be systematic- 
ally varied as well. 

The phase and gain analysis was performed by driving a unit 
magnitude sinusoidal signal at frequency R as the system output, 

y(t) = sin(fit) (3.15) 

To estimate the output, a second order LMS filter was implemented at 
sample rate <o s . The filter equations and parameter updates 
became 


y(k) = a i (k)y(k-l ) + a 2 (k)y(k-2) (3.16) 

e(k) = y(k) - y(k) (3.17) 

a l (k+1 ) = a x (k) + 2y£(k)y(k-l) (3.18) 

a 2 (k+D = <* 2 ( k ) + 2ye(k)y(k-2) (3.19) 


The system signal was initiated and the LMS filter transients 
were allowed to die for 20 seconds or 3 periods of system signal. 
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whichever was greater. The in-phase and quadrature components were 
computed, respectively, for one cycle as: 


P = — / y(t)sin(ftt)dt 

TT 

T 

0 = — / y(t)cos(ftt)dt 
7T 


(3.20) 


(3.21) 


The phase anrl gain were computed by the following relations 


<f> = tan 1 (^-) 


(3.22) 


M = 10 log in [0 ? ' + P 2 ] 


(3.23) 


The RMS error of prediction was also computed during the period of 
phase and gain evaluation. 

Figure (III-12) is a plot of phase and gain as a function of 
freciuency at a sample rate of 25 Hz. The phase maintains good 
agreement ( | $ | < 13 degree) up to an input signal frequency of about 
30 rad/sec. The gain remains relatively flat up through about 20 
rad/sec. Figure (TIT-13) shows a plot of the RMS error of predic- 
tion, confirming the region of difficulty at frequencies greater 
than 20 rad/ sec. 

As expected , it was observed that the sample rate is a prime 
factor in determining the frequency where the performance of the 
filter begins to degrade. Figure (TTT-14) is a plot of the ratio of 
frequency for a 15 degree phase margin to sample rate as a function 
of sample rate. A ratio of 1/4 to 1/10 hounds most of the answers, 
so it seems that good design practice would call for using a sample 
rate of 10 times the highest frequency being modeled by the LMS 
adaptive algorithm. 



PHASE AND GAIN PLOTS OF LMS FILTER 



FREQUENCY, ft, rad/sec 

Figure III-12 Phase and gain plots of second order LMS filter as 
a function of input signal frequency at a sample 
rate of 25 Hz. 


RMS ERROR OF PREDICTION FOR LMS FILTER 



FREQUENCY, ft, rad/sec 


Figure III— 13 RMS error of prediction for second order LMS 

as a function of input signal frequency at 
a sample rate of 25 Hz. 
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FREQUENCY RATIO FOR LESS THAN 15 DEGREE 
PHASE MARGIN AS A FUNCTION OF SAMPLE RATE 



SAMPLE FREQUENCY. Hz 


Figure III-14 — Frequency ratio for satisfying 15 degree phase 

margin for second order LMS filter as a function 
of sample rate. 
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B. MIMO CHARACTERISTICS OF IMS ALGORITHM 


A number of simulations were performed studying multi-input, 
multi-output identification to verify the multivariable formulation 
of the LMS algorithm developed in section II-C. No unusual charac- 
teristics were observed; as the number of parameters being identi- 
fied increased, the accuracy decreased as predicted. 

As previously mentioned, the mean parameter error increases 
linearly with the number of parameters being identified. There are 
2 parameters for the first order system of subsection III-A.l, 4 
parameters for the second order SISO system of subsection III-A.3 
and there are 10 parameters for a second order system with 2 inputs 
and 2 outputs. There are n+nrm parameters in a MIMO system of 

order n with m inputs and r outputs. Clearly, the number of 
parameters grows significantly for MIMO systems. 

During this research, systems were studied ranging from first 
order to tenth order; some were SISO and some had 2 inputs and 2 
outputs. Although the results are somewhat inconsistent (no effort 
was made to systematically match sampling rates, duration of adap- 
tion, plant dynamics or step size factors) a plot of average para- 
meter error is shown as figure (III-15). It is a plot of what is 
supposed as the best achievable for each problem. There is signifi- 
cant scatter in the data, but the anticipated trends are observ- 
able. One can expect a degradation in parameter estimation accuracy 
as the number of parameters increases. This is another good argu- 
ment for using the IMLMS formulation as it allows for fewer 
parameters to be estimated at each stage of the model building. 

Figure (III-15) also illustrates the fact that the numerator 
terras converge slower in time and have much larger parameter 
errors. Inaccuracy is substantially reduced by using the recursive 
LMS algorithm in conjunction with the batch least squares identifi- 
cation for the numerator terms as was developed in section II-E. 
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PARAMETER ERROR PERFORMANCE AS A FUNCTION OF THE NUMBER 

OF PARAMETERS 
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Figure III-15 — Average normalized parameter errors as a function 
of the number of model parameters. 


60 



C. CONVERGENCE IN PRESENCE OF UNMODELED MODES 


References [31-33,39-40,48-50] have indicated that LMS adaptive 
filters converge well even when used with mismatched model order. 
A simulation study to verify these assertions will be presented in 
this section. A fourth order SISO plant is modeled as equations 
(3.9) and (3.10) with 


A = 


A^Wj,^) 0 0 

0 0 0 0 
0 0 A 2 (t» 2 ,Z 2 ) 


b.[M 

L B 2 J 


H - [h l ,H 2 ] 


The submatrices are defined as 


Wi> = 


-“i " 2? i“i 


B i = 1 

L ij 


H i = [\ 0] 


(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 


Measurement and control influence terms are normalized by the first 
mode to define 


h = 



(3.29) 
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(3.30) 



The input signal is the pulsed wave of equation (3.11) where 
u amp “ *5, ft = 1 rad/ sec and the sample rate, oj s , is 20 Hz. 
The simulations were run for 200 seconds with p = .008 and for all 
runs wj = 1 rad/ sec and hj « bj_ = 1. 

The LMS adaptive filter was implemented assuming only one mode 
and the parameters of the second mode were varied. Identification 
accuracy was found by comparing the experimental LMS coefficients 
with the theoretical discrete coefficients of mode 1. A summary of 
the results with ?i = ^2 = *5 is plotted as figure (III-16). 

The total percentage error for the denominator terms is shown 
in figure (HI-16) as a function of the frequency of the unraodelled 
mode, 0 ) 2 * The LMS filter does remarkably well for 0)2 greater 
than 2. Looking at the curves with h = 2 and greater, it is appar- 
ent that the LMS filter tries to identify the lowest frequency mode 
until the signal power of the new mode becomes prohibitively 
strong. As the frequency of the unmodeled mode approaches the fre- 
quency of the modeled mode, 1 rad/ sec, the parameter error grows. 
When h > 2, the parameter estimate near 1 rad/ sec is nearly 

divergent as the LMS algorithm begins to follow the second mode. 

A survey was done varying and £2 but making them 

equal, and the same shape curves as figure (HI-16) were obtained. 
However, the actual parameter estimation accuracy was dependent upon 
the input sequence used to excite the system. With the same input 
as was used for generating figure ( III-16) and with varying the 
damping ratios between .005 and 1.0, the asymptotic total parameter 
estimation error varied between .12 percent and 5.9 percent. 

If the damping ratios of the two systems are not held constant, 
the impact of the input signal is felt more strongly. However, the 
same trends are apparent. If the spectral separation between two 
modes is wide and they have nearly equal signal power, the LMS 
filter can be expected to perform well in estimating the values of 
the parameters of the lowest frequency mode. 
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IMPACT OF UNMODELLED MODES UPON PARAMETER 
IDENTIFICATION ACCURACY 



Figure III-16 — Total parameter estimation error versus the 
frequency of the unmodeled mode. 



D. IMLMS CONVERGENCE CHARACTERISTICS 


The Incremental Mode LMS (IMLMS) algorithm was derived in 
section II-F and is useful because only part of the model is 
adapted. It allows the known part to be held constant while the 
unknown part is adjusted. This provides a natural way to perform 
model building — that is, by adding one mode at a time. Although no 
additional assumptions were needed to prove convergence of the IMLMS 
filter compared to the LMS algorithm, it is still of interest to 
verify the performance of the parameter estimates. 

The generic sixth order plant is modeled by assuming the 
linear, time-invariant, state space representation of equations 
(3.9) and (3.10), with A in modal form: 


A - 


h. 


B = 


L B 3 


H = [H^ H 2 , H 3 ] 


(3.31) 


(3.32) 


(3.33) 


The submatrices are defined by equations (3.26-3.28). For this 
simulation the plant parameters were selected as follows: 

^1 2 3 = 1*>^*>8* ras/sec 


b l ,2 ,3 h l,2,3 1 *» 2 * 5 > 1 


? 1,2,3 * 1 * * 1 * * 1 
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The sample rate, u> s , is 10 Hz so that <d s > 0)3 and the 
input signal is a pulsed wave given by equation (3.11) with 
u amp = * 5 and n=1 rad / sec. 

Time histories for the simulation are shown in figures (III-17) 
through (III-20). It is assumed that the first mode is known 
exactly and the IMLMS algorithm is used to identify the second 
mode. The third mode is not modelled and is treated as process 
noise similar to the previous section. Since the spectral separa- 
tion between modes 2 and 3 is 4 rad/ sec or a ratio of 1 to 2 , it is 
expected that good parameter estimation is possible. The output, y, 
the control input, u, and the modal positions, xj , X 3 , and X 5 
are plotted in figures (III-17), (III-18) and (III-19), 

respectively. 

The recursive parameter estimates for the two denominator 
parameters of the second mode are shown in figures (III-20) and 
(III-21). The parameters obtained 90 percent of their asymptotic 
limits within about 3 or 4 periods of oscillation of the second 
mode. The total parameter error was small (.23 percent). There was 
no tendency toward a bias and it converged quickly, even in the 
presence of another mode. So as predicted, this illustrates para- 
meter convergence properties similar to the LMS algorithm for the 
same number of parameters being adapted. 



TIME HISTORY OF OUTPUT 



TIME, t, sec 

Figure III-17 — Time history output of third order system used for 
studying IMLMS algorithm. 


CONTROL INPUT TIME HISTORIES 



TIME. t. sec 

Figure III-18 — Time history of control input used to excite third 
order system for studying IMLMS algorithm. 
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POSITION OF EACH OF 3 MODES 



TIME. t. sec 


Figure III-19 — Time history of the positions of each of three 
modes during simulation used to study IMLMS 
algorithm. 
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ADAPTION TIME HISTORY OFa 1 



Figure III-20 — Time history of denominator coefficient, 04 , for 
IMLMS algorithm demonstration. 


ADAPTION TIME HISTORY OF 0C 2 



Figure III-21 — Time history of denominator coefficient, <*2» for 
IMLMS algorithm demonstration. 
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Chapter 17 


STATE SPACE REALIZATIONS FROM MIMO TRANSFER FUNCTIONS 
A. PROBLEM DESCRIPTION 

The multi-input, multi-output (MIMO) version of the LMS and 
IMLMS adaptive filters yield the coefficients for the following 
discrete ARMA model in indicial notation 

A 

y^Oc) = a^Ck-i) + 0 u^Ot-i) (4.1) 

This represents a set of discrete transfer functions [21] 


H(z) 


11 4 

l B.z _i 
i=l 


1 


l 


i=l 


a 


-i 

i 


(4.2) 


The compensator portion of the adaptive control scheme proposed 
here will be implemented on-line using optimal control techniques. 
This requires a state space realization of the plant model to take 
advantage of the computer codes that are available for optimal 
design of linear quadratic gaussian (LQG) systems [51]. 

Finding a state space realization given a single-input, single- 
output (SISO) transfer function is straight forward [52]. However, 
finding a state-space realization given a MIMO transfer function as 
a state space representation is considerably more complicated. 
Several methods have been proposed [53,54], but they tend to be 
effective only when the coefficients of (4.2) are exact. The 
coefficients of (4.2) as determined by the LMS or IMLMS algorithms 
will not be exact. Additionally, a truncated, equivalent system of 
a higher order system may be desired. 
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A robust algorithm for finding minimal state space realizations 
was proposed in reference [ 55 ] • However, it is costly in terms of 

computational burden as it requires nonlinear programming tech- 
niques. In section C an improved version of the algorithm is devel- 
oped which greatly reduces the required computations. A further 
refinement is proposed in section D for the special case of a single 
mode with 2 inputs and 2 outputs (a 2x2x2 system). A partitioned 
linear algorithm is derived in section E which requires no itera- 
tions. In the final section of this chapter, some examples are 
given. 


B. ALGORITHM FOR STATE SPACE REALIZATIONS 

The algorithm of reference [55] is developed below in the 
s-domain, but it is easily extended to the z-domain. It converges 

only for stable systems (poles inside the unit circle). 

The transfer function matrix is represented as 


P(S ) = N ill 

* KS) d(s) ’ 


(4.3) 


and a state space realization is desired 


P(s) = H(sl-A) *B + L , (4.4) 

where A is a nxn matrix, B is a nxm matrix, H is a rxn 

matrix and L is a rxm matrix. P(s) is proper and rational and 
d(s) is the least common denominator of order n. 

L is found in the usual way 

L - p(s> . (4.5) 


The ARMA representations in this research are usually assumed to be 
of unit delay, so that L = 0. Define the following 
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M(s) = N(s) - Ld( s) 


(4.6) 


therefore. 


M( s) = H(sI-A) -1 Bd(s) (4.7) 


The denominator can be represented as 


d(s) = a rt s n + a s 11 * + a s n ^ + ... + a s^ (4.8) 

0 12 n 


or by the summation 


d(s) = l a s 1 (4,9) 

i=0 

The term (sI-A) - * can be represented by the exponential 
series given by [52] 

CO 

(si- A) -1 = l s _k A k (4.10) 

S k=0 

Substituting (4.9) and (4.10) into (4.7) yields 

oo Y1 

M(s) = “ H l s k /S l a .s* (4.11) 

5 k=0 i=0 n 1 

which also equals 


rj oo 

M(s) - H l l s 1-k-1 A k a B (4.12) 

i=0 k=0 n " 
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If the following definition is made 


k’ = -k+i-1 , (A. 13) 

then 

k = -k'+i-l , (A. 14) 


and 


M(s) = H l l 

i=0 -k'+i-l=0 


k' -k+i-1 
s A a 


n-i 


(A. 15) 


Adjusting the k' summation gives 


M(s) = H l 

i=0 k'=i-l 


r k' -k'+i-l 
), s A a j B . 


n-i 


(A. 16) 


Breaking (A. 16) into positive and negative summations of k' gives 


n 0 


M(s) = H l l s A 

i=0 k'=i-l 


k' -k'+i-l 


n-i 


? “r k' -k'+i-l 
+ l I s A a , B 

i=0 k'=-l 


n-i 


(A. 17) 


Note that a zero i gives a negative k' . So in the positive k’ 
summation, i goes from 1 to n; and in the negative k' summation, 
i goes from 0 to n. The negative summation of equation (A. 17) 
equals 

1 s^’- 1 l A i a n _ i (A. 18) 

k’=-l i-0 
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But 


l 


1=0 



la + Aa + 
n n-1 


A n 

A a r 


(4.19) 


Equation (4.19) Is equal to zero by the Cayley-Hamilton Theorem 
[52] since d(s) is the minimal polynomial of P(s). Therefore, 
equation (4.17) becomes 


n i-1 


M(s) = H l l s A' 
i=l k'=0 


k' -k'+i-l 


a . B 
n-i 


(4.20) 


Equation (4.20) is solved iteratively by specifying a set of 
residuals, R(s), for minimization using nonlinear programming 
techniques. 


n-1 

R(s) = M(s) - l 
k=0 


k V n-l-i 

s I a .HA B 
i=k 


(4.21) 


Ideally the elements of R(s) approach zero at a solution determined 
by the numerical procedures outlined below. 

Developing the following indicial notation 


R(s) = R n s n_1 + R.s n “ 2 + ... + R ,s° 
0 1 n-1 


(4.22) 


r T • • • T • • • T 

ill il2 ilj ilj 



r t* • • • r • • • r 

irl ir2 ikj irm 


(4.23) 


and if M utilizes the same notation, we can equate indices in 
(4.21) and write 
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(4.24) 


n-i n-i n-i 

8 r ikj = 8 m ikj " 8 


n-1 


r tl . n-l _ £ 

i a, .HA B 

£=1 X ' i 


It is desired to make the sum of the squares of the residuals, 
r^kj, as small as possible. Since values of m ikj ma y vary 

significantly in magnitude, causing numerical problems, (4.24) is 
normalized by m ikj (in the computer code, m ikj ^- s se t to 1 if 
less that 10“^). Normalizing and equating terms of like powers of 
s n ~* gives 


"ikj - X 


ikj 


m 


ikj 


(4.25) 


The implied summation limits for i, j and k are 0 to n-1, 

1 to m, and 1 to r. Hence, (4.25) implies that nrm equations are 

required for the minimization of r^^j, 

A unique solution of (4.20) is obtained through the minimiza- 
tion of r ^ikj from* (4.25) by choosing an A matrix with the 
proper eigenvalues from the coefficients (4.9) and by specifying one 
nonzero element in each row of B or each column of H. The 

remaining elements of B and H are found iteratively using non- 
linear programming techniques. 

For convenience, the A matrix is chosen in observer canonical 
form [21,52] and the elements of the first row of H are chosen as 
unity. Equation (4.22) is solved computationally by using Stewart’s 
adaption of the Davidon-Fletcher-Powell (SDFP) algorithm [56,57], 
where the gradient is computed by numerical differentiation. The 
minimization is performed over (mfr-l)n optimization variables 

using a computer code for the SFDP algorithm available at Langley 
Research Center. 

Initial implementations of this multivariable state space 
algorithm were found to be extremely sensitive to initial condi- 
tions. Two solutions were adopted to improve operating 
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characteristics under varying initial conditions: 1) If a satisfac- 

tory convergence is not obtained in a certain number of function 
calls, the routine restarts itself with a new set of parameters 
chosen randomly between some preselected boundaries; and, 2) an 
optimization variable transformation was used to scale the variables 
and apply constraints to the region that the search is allowed to 
take place over. 

The optimization variable transformation [58] from the para- 
meters, pj_, of H and B to the optimization variables, z^, is 
given by 


(U.-L.) _ U.+L. 

P, = o sin [— z.J + [ — - — J 


(4.26) 


where z-^ is the ith variable the SDFP code optimizes, and 

are the upper and lower values, respectively, that are con- 
straining the ith parameter, p£. The inverse transformation is 
clearly given by 


rU .+L . 'i 

P, ~ l i i J 

_ 2 . -If 2 'i 

Z i v Sin (IL-I^) J (4.27) 

2 

The key feature of this transformation is that the optimization 
codes search parameters over a range of -1 to +1 . Furthermore, any 
value of will always return a value for p^ in the range 

Li to Ui. Hence, the code optimizes variables of approxi- 
mately the same size and the Pi's have been effectively 

constrained to lie between Li and Ui* 

In the next three sections, further improvements in the 

algorithm are developed. Section F of this chapter gives examples 
illustrating the use of these algorithms. 
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c. MODIFIED ALGORITHM BY PARTITIONING 


An algorithm for obtaining state space realizations from 
multivariable transfer functions [55] was developed in the previous 
section. It is noted, however, that the solution of equation (4.20) 
through the normalized error residuals of (4.25) involves the use of 
nonlinear programming with the optimization code searching over 
(m+r— l)n variables. In the simplest useful case of a 2x2x2 system 
(nxrxm), there are 6 degrees-of-freedom the program must search 
over. A typical representation for the longitudinal motions of an 
airplane would require 12 parameters (4x2x2). The overall adaptive 
control algorithm should be able to translate the parameters 
estimated by the LMS or IMLMS filters into state space realizations 
in real-time on microprocessors prior to the suboptimal compensator 
design. Performing a search over this many parameters would be 
quite slow. 

In this section an improvement to the algorithm of [54] is 
proposed that reduces the execution time by reducing the number of 
degrees of freedom. The approach used is to partition the problem 
into a linear part and a nonlinear part. The linear part is solved 
with simple matrix operations using the current values of the 
iterated nonlinear part. In this fashion the size of the nonlinear 
portion of the problem is reduced. 

Take the matrix equation of (4.25) for the s° term of the 

numerator (i=0) 


n 


B 0 


M - I a HA 
0 1=0 


n-1 — l 


B 


(4.28) 


Since it is desired that Rq approach zero, we can rewrite 
(4.28) as 


- l a^-^B 

0 Jt=0 


(4.29) 
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which is just a different form for one power of s from (4.20) 
Assume H is known, then (4.29) can be solved for B 


B = 


I 


1=0 




n-l-£ 



(4.30) 


The () + indicates a matrix pseudo inverse, which is needed if r 
and m are not equal. In the likely event that r and m are 

equal, then a simple matrix inverse is needed. 

Equation (4.30) implies that once the elements of H have been 

chosen, B can be computed using linear algebra. Therefore, the 

nonlinear optimization code now only has to search over the elements 
of H; and, B is computed by (4.30). Hence, the solution of 
(4.20) has been effectively partitioned into a nonlinear part (find- 
ing the remaining elements of H) and a linear part (finding the 

elements of B) . 

The degrees of freedom for the algorithm of [54] have been 
reduced from (r+m-l)n to (r-l)n. In the case of a 2x2x2 
example, this corresponds to a reduction from 6 degrees-of-f reedom 
to 2, which would translate into more than a 66 percent savings in 
execution time for the optimization code. The 4x2x2 example is 
reduced from 12 to 4 parameters. It is anticipated that this would 
correspond to a reduction of nearly 85 percent in execution time. 
These savings are significant, especially when on-line implementa- 
tions are considered. Additionally, the accuracy of the algorithm 
can be expected to improve with a reduction in the number of free 
parameters. 

It may be advantageous to alternatively specify n nonzero 
elements in B and find the remaining elements in B by nonlinear 
optimization. Equation (4.29) can then be rearranged for a linear 
solution of H 

H « M L ( l a A n_1 “ Z B) + (4.31) 

£=0 
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Using (4.31) would result in (m-l)n degrees-of-freedom. In 
this research (4.30) was used, but the choice should be based upon 
the relative magnitudes of ra and r. 

D. PARTITIONED ALGORITHM WITH ANALYTICAL GRADIENTS 

FOR SINGLE MODE 

In the previous section, the algorithm of [54] was partitioned 
into a linear and a nonlinear part to improve computational effi- 
ciency. The implementation of the algorithm with Stewart's adaption 
of the Davidon-Fletcher-Powell (SDFP) computer code utilizes 
numerical differencing techniques for the computation of the gradi- 
ents. This is time consuming and is required only when analytical 
relations for the gradients are not available. General relations 
for the gradients were not possible because of the inverse/pseudo 
inverse aspect of equation (4.30) that is used in the quadratic form 
of (4.25). However, it is possible to derive analytical gradients 
for special cases. 

In the multivariable examples used in this research, the prob- 
lem of adding a single mode with two inputs and two outputs (2x2x2) 
to a model was encountered frequently. In order to speed up the 
partitioned algorithm further, analytical relations for the first 
partials and second partials of the cost function with respect to 
the parameters of H were found. This allowed direct implementa- 
tion into Newton's method of optimization [22], The iteration in 
terms of optimization variables, z^, would be processed for the 
kth iteration as 


*i< k+1 > - z i< k > + j 22 ‘ lj 2 


(4.32) 


where J is the cost function defined from (4.25) as 


J 


n-1 m 

l l 


i=0 j=l 


l 


k=l 



(4.33) 
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j z and J zz are the first and second partials of J with 
respect to the optimization variables, z-^. 

For the 2x2x2 case considered here, B is found from (4.30) as 

B = (a 1 H+HA) -1 M Q (A. 34) 


and the residuals come from (4.25) as only the s term since (4.34) 
exactly satisfies the i=0 term 


_ [m. ... -a-.HBV 

-2 L ljk 0 J 

r = 

ijk 2 

m. 
ijk 


(4.35) 


For convenience H is taken as 


H = 


1 1 


h 21 h 22. 


(4.36) 


so that the parameters that are sought in the nonlinear part become 


£ = 


21 


L h 22. 


(4.37) 


The first and second partials of J with respect to £ were 
computed on MACSYMA [59,60]. The 4 second partials and the 2 first 
partials took up over 50 lines of Fortran code. Although the compu- 
tational burden is lower, the complexity is much greater. 

The optimization variable transformation of [58] was used in 
this case, as well, to scale the parameters and to apply limits to 
the range the search is taken over. Hence, the partials of J with 
respect to z^ are computed as follows 
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(4.38) 


3j_ _ 9J_ . 9p l 
3z ± 3p ± 3z ± 


3 2 J _ 3 2 J # 9p i 

3z.3z. 3p.3p. 3z. 

i J i J i 

The correction terms in (4.38) and (4.39) for transforming the 
gradients between the z-domain and the p-domain are developed from 
(4.26) as 


!z j 

3z . 


(4.39) 


8p i (U i _L i ) TT fit , 

3z ± 2 2 COS 4 V 


(4.40) 


This 2x2x2 version of the multivariable state space realization 
had significant savings as the number of required function calls was 
reduced 50 percent over the partitioned method of the previous 
section. 

In this section the analytical gradients, Jp and Jpp> 

were found on MACSYMA [59,60] only for the 2x2x2 problem. MACSYMA 

was used to directly obtain Fortran code, indicating the possibility 
of finding the gradients for any specific problem if a particular 
application warranted it. For example, a 4x2x2 system would be 

useful for airplane applications. So the general procedures 
developed in this section would be useful for developing adaptive 
control systems for other cases. However, the lack of generality 
and the length of the sensitivity equations limit the widespread 
applicability of this approach. 

E. PARTITIONED LINEAR ALGORITHM 

After performing the simulations for this report, another 

method was found which partitions the problem into two linear 
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problems, 1 One is solved with a simple matrix inverse and the other 
is solved with a pseudo inverse (least squares). No nonlinear 
programming techniques are required. This approach is derived here 
and results are presented in the next section which illustrate the 
increased efficienty and accuracy. 

The partitioned linear algorithm is derived by rewriting the 
infinite matrix polynomial equation (4.10) as a finite partial 
fraction expansion 


( sI-A) 


I 5 C n- 
£=0 
d( s) 


n-£ 


(4.41) 


where C n _£ are the numerator matrices. Then by equating the 
transfer function numerator matrices in powers of s, equation 
(4.41) can be partitioned into two linear problems. 

The numerator terms from (4.6), (4.7) and (4.41) are equated as 

n- & * n 1 « 

H l s C B = l s M . (A. 42) 

n r\ n Ao n r\ n Af 


As for the previous approaches, a nonzero element is specified in 

each column of H. For convenience and compact notation, l’s are 

chosen for all the elements of the first row. If the notation 

(,)^ is used to indicate taking the first (uppermost) row of (.) 
and (.) L is used to indicate taking the remaining (lower) rows 

of (.), H can be partitioned as 


H = 



(4.43) 


^his method was proposed to the author by Professor Arthur E. 
Bryson, Jr., of Stanford University in a personal letter dated April 
20, 1983. This was later refined in a letter from Prof. Bryson on 
May 14, 1983. The author has subsequently generalized the approach 
for high order systems for any number of inputs and outputs. 
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where 


H U = [1 1 ... 1] (4.44) 

with dimension lxn. has dimension (r-l)xn. The desired 

numerator matrices, M(s), are partitioned in the same fashion 


M = 



(4.45) 


The dimension of is lxm and the dimension of is 

( r-1 )xm. 

Define a new matrix, T H > which stacks the first rows of the 
matrix products 


T = 
H 


r H U C 

hV 


‘-H U c 3 - 


(4.46) 


Also, the first rows of the M(s) matrices are stacked vertically in 

V B 


r<v u i 

(m 2 ) u 

_(m 3 ) u _ 


(4.47) 


Equation (4.42) is partitioned using (4.46) and (4.47) into a 
linear equation in terms of B. 

T r B = V fi (4.48) 
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B is found by premultiplying (4.48) by Tjj 1 


B = 



(4.49) 


The remaining nr-n values in H are found by solving the 
remaining equations of (4.42). Define two new matrices as 


T b = [CjB C 2 B ... C n B] 

(4.50) 

■ [(M x ) L (^ ) L ... (M n ) L ] . 

(4.51) 


The dimensions of Tg are nxnm and the dimensions of Vg 
are (r-l)xnm. The lower partition of (4.42) is written as 

H L T b = V R , (4.52) 


which is a set of nmr-nm equations for nr-n unknowns. The unknown 
values for H are found in a least squares sense as 

h l = VV + (4.53) 

where () + indicates a pseudo inverse. 

The numerical solution is found in a direct way with no 
requirement for nonlinear programming. Nonzero values for the first 
row of the output distribution matrix, H, are assumed. The partial 
fraction expansion for the numerator terms of (sI-A) - l, C, is com- 
puted using the Leverrier-Souriau-Faddeeva-Frame formula [51 ,52] . 
The control influence matrix, B, is computed using (4.49). Then a 
least squares solution for Hl, the remaining terms of the output 
distribution matrix are computed using (4.53). This approach 
requires no iteration and is computationally fast. The most strin- 
gent computational burden is the need for two matrix inverses of 
size nxn. 
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F. NUMERICAL EXAMPLES 


In this section a 2x2x2 transfer function from [54] is 
considered and the algorithms of the previous four sections are used 
to find a state space realization. The answers are compared by 
evaluating the residual at the solution. First a problem with exact 
coefficients is given and secondly the same problem with inexact 
coefficients is utilized. The exact transfer function is 



The chosen A matrix becomes 



and the M matrices are 



The H matrix is assumed to be 



(4.54) 


(4.55) 


(4.56) 


(4.57) 


(4.58) 


Initial guesses for the parameters of the partitioned method were 
h2i=l and h22 = *25. 

The exact solution for the state space realization is 


B 



(4.59) 
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A comparison of the straight SDFP, the partitioned SDFP, the 
Newton 1 s iterative, and the partitioned linear methods are shown in 
table (4.1) . The partitioned SDFP algorithm reduced the number of 
function calls from 25 to 7 and improved the accuracy of the fit. 
The partitioned Newton’s method with analytical gradients reduced 
the number of function calls further to 3. The partitioned linear 
method required no iterations and was much more accurate than the 
other three. 

Table 4.1 — Convergence characteristics of realization algorithms 
for a problem with exact coefficients and of size 
2X2X2. 


METHOD 

NO. OF FUNCTION CALLS 

RESIDUAL 

SDFP 

25 

4.01 E -07 

PARTITIONED SDFP 

7 

1.21 E -14 

PARTITIONED NEWTON 
WITH ANALYTICAL 
GRADIENTS 

3 

1.70 E -11 

PARTITIONED 

LINEAR 

LEAST SQUARES 

1 

2.31 E -18 

















The solution that the partitioned linear method achieved was 


B = 


-2.39x10 


-10 


1.0 


.25 

.75 


(4.61) 


H 


1.0 1.0 

4.0 9.57x10 


(4.62) 


Now consider the case where the transfer function is an inexact 
representation of the same system. This might occur with experi- 
mental data, for example. The inexact transfer function is given by 


P(s) 


”l .lls+6.19 1.028+3.25“ 

_ 4.11 . 9 s+3 . 1 2 _ 

s 2 +5.03s+6.62 


(4.63) 


The matrices A, M| , M 2 and H are defined in the same 
fashion as equations (4.55-4.58)* The results of applying the three 
algorithms is shown in table 4.2. Partitioning the SDFP algorithm 
resulted in a reduction of function calls from 39 to 10. Computing 
analytical partials of the cost function allowed a further reduction 
in function calls to 5. The partitioned linear method gave accept- 
able performance with no iterations. 

The solution that the partitioned linear method obtained for 
the inexact case is given by 



(4.64) 


(4.65) 
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Table 4.2 — Convergence characteristics of realization algorithms 
for a problem with inexact coefficients and of size 
2X2X2. 


METHOD 

NO. OF FUNCTION CALLS 

RESIDUAL 

SDFP 

39 

3.40 E -07 

PARTITIONED SDFP 

10 

4.44 E -06 

PARTITIONED NEWTON 
WITH ANALYTICAL 
GRADIENTS 

5 

7.36 E -07 

PARTITIONED 

LINEAR 

LEAST SQUARES 

1 

3.51 E -4 


Using (4.64) and (4.65), P(s) is 


l.lls+6.19 1.02S+3.25 

♦057+4.086 .907s+3.15 

s 2 +5.03s+6.62 

Most approaches to state space realizations of this second 
example would have resulted in a fourth order system due to the 
inexactness of the experimentally-determined coefficients. Since 
the LMS or IMLMS algorithm is implemented by stating the desired 
size of the block to be added, it is advantageous for the algorithm 
to return a minimal state space realization in a least squares, best 
fit sense. If the value of the residual at the solution is too 
large, further estimation of the parameters by the output error 
identification scheme would be warranted. 



6 _ 


(4.66) 
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Since real-time applications for the algorithm using micropro- 
cessors are anticipated, the linear partitioned method appears to be 
the most promising. It requires no iteration and is totally general 
for any size system with any number of inputs and outputs. Addi- 
tionally, it is coraputationlly simple as it requires only two matrix 
inverses of matrices with rank equal to the system order. It sacri- 
fices some accuracy for inexact problems, but these differences were 
negligible for the cases considered. 
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Chapter V 


SUBOPTIMAL COMPENSATOR DESIGN AND IMPLEMENTATION 
A. INTRODUCTION 

Section I-C describes the approach to adaptive control in this 
research. A fast, robust identification scheme estimates a set of 
discrete transfer functions for the process being controlled (chap- 
ters II and III). These transfer functions are then converted into 
a state space realization (chapter IV). Using the compact matrix 
notation of the state space realization, digital compensator logic 
is designed using linear-quadratic-gaussian (LQG) techniques (this 
chapter) . 

Figure (V-l) is a block diagram that shows the form of the 
suboptimal compensator. Albeit LQG or so-called optimal control 
techniques are used, the compensator is referred to as suboptimal 
because generally no prior knowledge of noise covariances is avail- 
able and because the model utilized for the compensator design is 
understood to be an inexact representation of the plant. Constant 
gain matrices K, C, and F are the output of the design pro- 
cess. To keep the computational burden low for the digital compen- 
sator, both the filter and the controller are implemented with their 
steady state values. The LQG digital designs were performed using 
ORACLS [51]. A more complete description of the computing tech- 
niques is given in appendix I. 

Some simplifications were required to obtain the capability to 
perform on-line designs. The form of the weighting matrices was 
assumed; and, it was assumed that the separation theorem [52] holds 
in spite of plant-model mismatch. This allows the filter and the 
controller to be designed sequentially. However, the overall 
system is checked for stability of the compensator. 



SUBOPTIMAL COMPENSATOR BLOCK DIAGRAM 



DIGITAL COMPENSATOR 


Figure V-l — Block diagram showing fundamental structure of 
suboptimal discrete compensator. 
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B. STEADY STATE KALMAN FILTER DESIGN 


Once the order and parameters of the model have been found, an 
asymptotic Kalman-Bucy filter [63] is found using ORACLS [51]. The 
system is modeled as a discrete, LQG, time- invariant process. The 
state equation is given by 

x(k+l) = *x(k) + T w(k) (5.1) 

— W ■ 


with output 


^(k+1) = Hx(k) + v(k) (5.2) 

where $> is the state transition matrix > r w is the discrete 
influence matrix of process noise input, w. It is assumed for 
simplicity in most filter designs that F w , is an identity matrix 
of order n* The elements of $ and H are determined by the state 
space realization methods of the previous chapter from the identi- 
fied parameters of the LMS or IMLMS algorithms* 

The process noise w(k) and the measurement noise _v(k) are 
modeled as purely random gaussian vectors with zero means and co- 
variance matricies Q and R. The solution requires that w(k) and 
v(k) are mutually uncorrelated. Since prior knowledge of the pro- 
cess is limited due to the model building approach being proposed, 
the on-line design is implemented by assumming that Q and R are 
identity matrices of rank n and r, respectively. The designer 
choices p, a scaler multiplicative factor of Q, as a means of 
varying the relative weights in the cost function. This simplifica- 
tion is possible since the output variables are scaled during the 
identification process. 

The cost function is defined as 

J = lim e[e^ ( k)We(k) ] (5.3) 
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where _£ is defined as 

e(k) = _x(k) - jt(k) (5.4) 

and W is a weighting matrix which does not appear in the 
computation. If the pair ($,r w ) is stabilizable and the pair 
($ ,H) is detectable, then a solution to the optimal observer problem 
exists and the predicted state, is given by 


x(k+l) = ^(k) + K[x(k)”Hx(k)] , (5.5) 


where 


x(0) = E[x(0) ] . 


(5.6) 


The steady state filter gain is found from 


K = $PH T (R+HPH T ) 1 , 


(4.7) 


where P satisfies the discrete algebraic Riccatti equation 


P = ($-KH)P('I>-KH) T + KRK T + PQ . (4.8) 

The matrix P represents the constant, steady-state variance 
matrix of the reconstruction error, e(i), given by 

lim E[x(k)x T (k)] = P (5.9) 

i=0+-°° 


The product of this computation is K, the Kalman filter gain 
matrix in figure (V-l). The designer chooses the weighting factor, 
p, and the on-line code uses the values of $ and H to find a K 
matrix. Some skill is needed in choosing p and it could 
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conceivably be delegated to the heuristic, problem solving logic of 
the adaptive processor. Clearly, if this approach was being imple- 
mented on a process where prior knowledge can be used to estimate Q 
and R, an optimal Kalman-Bucy filter may be possible. 

C. STEADY STATE REGULATOR DESIGN 

After the order and parameters of the model have been deter- 
mined, a state space representation found, and an asymptotic filter 
gain computed, an asymptotic quadratic regulator gain, C, is calcu- 
lated (see fig. (V-l)). The discrete, time-invariant, linear, 
optimal output regulator problem is solved using ORACLS [51]. The 
system is given by 


x(k+l) = $jx(k) + TuCk) » (5.10) 

and the output vector is 


jK k) = Hjc( k) . (5.11) 

The cost function that is optimized is 

N “ X . T T 

J = lim E l [^(k+OpQ v(k+l) + u (k)Q u(k)] (5.12) 

N-k» i=0 “ 

subject to the constraints of (5.10) and (5.11). The matrices 4>, 
T and H are obtained from the estimation and state space 
realization sections of the code. Qy and Q u are weighting 

matrices that are assumed to be identity matrices of order r and 
m, respectively. The design parameter, p, is a scalar multiplica- 
tive factor of Qy 

If the pair ($,T) is stabilizable and the pair ($,H) is 
detectable, a solution to the optimal control exists and is given by 

u(k) = -Cx(k) , (5.12) 
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where 


c = (R+r T pr) -1 r T p$ 


(5.13) 


and 


P = ($-rc) T P($-rc) + C T Q u C T + H T pQ y H (5.14) 

The steady-state controller is designed by computing the 

constant regulator gain C. The code uses as input the matrices 
T and H and the designer chooses p. The weighting factor, p, 
is chosen by the designer as a relative trade in the cost function 
between output power and allowable control input power that is mini- 
mized. The convenience of utilizing p as the design parameter is 
possible since the variables are scaled during the identification 

process. It may be desirable to allow the adaptive processor the 
opportunity to heuristically select the magnitude of p. 

The separation theorem was used to design the controller and 

observer separately, but it should be noted that when actually 

implemented they are linked by 

^i(k) = -Cx(k) (5.14) 

That is, the estimated states are fedback. Stability is not guaran- 
teed since the transition matrices of the estimator and regulator 
may not be sufficiently accurate representations of the process 
being controlled. This assumption that the separation theorem holds 
is not necessarily true, but the suboptimal compensators designed in 
this fashion do exhibit good control characteristics for a wide 
range of problems. 
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D . NON-ZERO SET POINT FOR DISCRETE REGULATOR 


During the model building process proposed here, it is 
necessary to sufficiently excite all the modes of the system. 
It was observed in chapter III that square wave pulses of the 
control inputs are a good way to do this. Although it may not be 
desirable during a control task to indiscriminately pulse the 
controls, some sort of dither signal is necessary to insure the 
convergence of the identification schemes. 

An alternative approach, which has proven to he successful, is 
to command the output to follow a pulsed square wave pattern. If 
the plant model is good, the output will closely follow the. com- 
manded output. However, if there are unmodeled modes present, thev 
will most likely he excited and discerned in the outnut. The reason 
for this is that feedforward gain matrix, F, is based on the plant 
model. 

The feedforward gain matrix, F, for a discrete regulator is 
based upon the non-zero set point calculations from reference (64]. 
Since the calculation is fundamentally requiring an asymptotic 
response to step commands, the zeros have been neglected [21], The 
state equations are 


x(k+l) = $x(k) + Iu(k) 


y(k) = Hy(k) + T,u( k) 


T he steady state equations occur when 


x = $x + Tu 
— ss — ss — ss 


v = Hx + T,u 

— SR — SS — SS 


Combining terms in (5.17) yields 


0 — ($-I)x + Tu 

— — ss — ss 


(1.15) 

(5.16) 


(5.17) 

(5. 18) 


( 5 . 19 ) 
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writing (5.19) and (5.18) in matrix form 



Solving for _Xss anc * Uss yields 


X 

— ss 


"*-I 

r" 

-1 

0 1 

u 

ss_ 


H 

L 


ksj 


(5.20) 


(5.21) 


The matrix inverse is partitioned as follows: 


$-1 r 

H L 

where the dimensions of S^ are nxn, S 2 are raxn, are 

nxr and F 2 are mxr. The steady state values of the states and 

control inputs are 


— l 

- Sl 

F. 

_ 

1 

1 


Csl 

CO 

1 

F 2. 


(5.22) 


x = F, 
— ss 1 


-ss 


(5.23) 


^ss = F 2 X 


ss 


(5.24) 


Since it is desired to have the steady state output approach 
the steady state commanded output, let 

Xss ’ Xc (5 - 25) 

Also, in the steady state x Qg -» -x g g> so from (4.14) 

u - -Cx . (5.26) 
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Adding a steady state component to both sides of (5.26) implies 
that corrective inputs would occur only for deviations from the non- 
zero set point 


Hss “ “ C( ^ss ) < 5 * 27 > 

Solving for the total control _u 

u = u - Cx + Cx (5.28) 

— — ss — — ss 

Substituting (5.23) and (5.24) into (5.28) gives 

u. “ F 2c “ Cx (5.29) 

where 

F = (CF 1 +F 2 ) (5.30) 

The constant feedforward gain matrix, F, is implemented in the 
suboptimal compensator of figure (V-l). The commanded output, 
y CO m» * s use d to pulse the output so as to excite the unmodeled 
dynamics. The amplitude and frequency of the square wave pulse can 
be selected by the designer. Alternatively, these design parameters 
could be systematically varied and optimized by the heuristic learn- 
ing logic of the adaptive processor. 

E. COMPENSATOR IMPLEMENTATION STRATEGIES 

The overall adaptive control scheme described in section I-C 
requires that the compensator design be updated periodically. The 
output error identification scheme is in parallel to the plant and 
the compensator (see fig. (V-2)). Information from the identifica- 
tion block passes to the compensator block only when the compensator 
is to be redesigned, which may include a change in model order as 
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OVERALL SCHEME FOR ON-LINE EQUIVALENT SYSTEM 
IDENTIFICATION FOR ADAPTIVE CONTROL 



IDENTIFICATION INPUTS 


Figure V 2 Block diagram showing on-line equivalent system 
identification scheme for adaptive control. 










well as model parameters. Information is passed when divergence of 
the parameter estimation scheme is detected, so divergence problems 
like the ones described in subsection I-B.2 do not occur for this 
approach. However, advantage is gained at the expense of adaption 
speed; it takes time to collect data, build a model, check for 
consistency and compatibility, and then to design a compensator. 

The questions addressed in this section are related to the 
issue of how often to update the compensator design. The strategies 
for choosing and evaluating figures-of-merit are also discussed. 
The fundamental issue is how long should the recursive LMS or IMLMS 
algorithms be run before a satisfactory convergence is obtained; 
and, once obtained, when is it warranted to replace the compensa- 
tor? Therefore, the scheme must have a means of self-evaluation so 
overall system performance can be enhanced. 

E— 1 Restructurable Control 

The first case considered is when the adaptive control scheme 
is used in a conventional way. That is, the order of the model is 
constant, but parameters have a significant uncertainty, time depen- 
dence or possible failure modes that might require the compensator 
to be redesigned (restructured) to augment the control system 
performance. 

Three primary figures-of-merit are used for evaluating subsys- 
tem performance. The parameter convergence of the LMS (or IMLMS) 
algorithm is initially determined by considering an RMS error of 
prediction from the recursive identification scheme, 


N 




1/2 


LMS 


N 


(5.31) 


where yi.Ms(k) is the predicted output from the LMS or IMLMS 
algorithms, N is the number of samples. 

If a LMS is a ^ ove some minimum threshold, say 10”^ in 

normalized units, it is assumed that convergence has not been 
achieved. If <?LMS *- s con *puted over a sufficiently long time and 
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is above some maximum threshold value, it is assumed that the LMS 
algorithm is diverging and should be reinitialized. 

The estimation subdivision of the compensator is evaluated 
using an expression similar to (5.31) 



where £(k) = Hx(k) is the predicted output of the Kalman filter 
portion of the compensator. 

The RMS estimation error, is a measure of how well the 

model utilized in the compensator design is doing in predicting the 
response of the plant. If Oj, after a certain minimum period of 
time, is above some preselected value, a new model with the para- 
meters obtained from the LMS or IMLMS algorithms is implemented. 
a-[ includes modeling errors, process noise and measurement 
noise. The threshold value for that is selected needs to be 

considered by the designer in the context of these other distur- 
bances. 

Another performance index which is critical for the restructur- 
able control problem is the RMS error of commanded output. It is 
computed by comparing the actual output with the commanded (desired) 
output. Typically, a set of square wave pulses is commanded. This 
f igure-of-merit is computed by 


a 

c 


k=l 

N 


(5.34) 


If a c becomes too large, either the model should be updated 
or the regulator portion of the compensator should be redesigned 
with a higher gain factor, p. The RMS error of commanded output, 
a c , includes the effort of the process noise and measurement 
noise, as well as the performance of the compensator. 
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Determining if a failure has occurred is a complicated and 
unresolved problem [65,66] that will not be considered here. How- 
ever, if a failure has been detected, the adaptive control logic 
proposed here for a restructurable control problem is applicable as 
shown in the flow chart of figure (V-3). This chart shows an 
approach that was successful for the examples considered; it con- 
sists of heuristic rules selected by the designer to implement an 
adaptive controller* The delay in time at the second decision level 
and the other checks prevent the recursively identified parameters 
from being fed directly through to the compensator, avoiding some of 
the divergences mentioned in chapter I. 

The logic of figure (V-3) is meant to systematically isolate 
what needs to be improved and is meant to operate simultaneously 
with the control, estimation and identification processes* If a 
failure is detected, a predetermined amount of time is allowed to 
pass so the identification schemes can have time to converge* If 
the control index is not satisfied, the model prediction is 
checked. If the model prediction is satisfactory, the controller is 
redesigned; otherwise, the convergence of the LMS algorithm is 
checked. Once a LMS ^ ow enou gh, the model parameters are 

updated, a state space realization is found and a compensator is 
designed using the techniques of this chapter. Once a compensator 
has been brought on-line, it may be necessary to further update the 
model or redesign the compensator as the LMS algorithm convergence 
improves. 

The prime limitation of this approach is that an unstable 
system could fail or break before the adaptive process could work. 
In general, the scheme has inherent delays to prevent adaptation 
divergence, but these may result in unsatisfactory performance. The 
next chapter shows an example of where the approach in this subsec- 
tion was used successfully to restructure a control system after a 
component failure. 

E-2 Model Building 

A prime contribution of this research is the IMLMS algorithm 
which includes the unstructured model uncertainty in the adaptation 
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DECISION LOGIC FOR RESTRUCTURABLE CONTROL 


START 



V-3 Flow chart showing decisions required for 

restructurable control applications of adaptive 
control algorithm. 


102 








process. Part of the model is assumed to be known and the other 
part is adapted by IMLMS algorithm. When satisfactory convergence 
of the IMLMS filter has been achieved, the order of the compensator 
model is increased as an additional mode is added. This process can 
be repeated over and over, thereby building a model. The key ques- 
tion is when to add a mode and when to stop. 

The model building logic is shown in figure (V-4); it uses the 
same f igures-of-raerit as the restructurable control application, 
°LMS> a I> and °L* However, in addition, the designer selects 
a minimum model order to begin control and maximum model order for 
the model building process. These parameters are chosen based on 
prior knowledge of the plant and its environment. Stability prob- 
lems may occur if a plant is controlled with too few modes repre- 
sented in the compensator design model. The maximum limit is a 
practical one which may be dictated by hardware constraints. 

If the maximum model order limit has not been exceeded, the 
adaptive model building logic tests to see if enough time has 
elapsed since the last control system modification, this allows the 
IMLMS algorithm enough time to converge. The designer chooses this 
based upon prior knowledge of system dynamics and the chosen sample 
rate. If a c is satisfactory, no adaptation is required. Other- 
wise, Oj is tested. 

If Oj is significantly bigger than the last time it was 
computed, experience has shown that this is because the model order 
was too large the last time a mode was added. Most likely the 
last mode was added based upon trying to identify nearly white 
noise, and when implemented in the Kalman-Bucy filter, degrades its 
ability to predict the output. A successful solution to this pre- 
dicament is to reduce the internal structure order by removing the 
last mode. 

The controller portion of the compensator is redesigned if a c 
is too large and is satisfactorily small. If is too 

large, another mode is desired for the compensator model. It is 
added once the IMLMS algorithm satisfies its convergence check. The 
controller portion of the model is implemented only if the number of 
modes in the compensator model is above the minimum value required 
by the designer. 
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MODEL BUILDING DECISIONS FOR ADAPTIVE CONTROL 


START 



Figure V-4 — Flow chart showing decisions required for model 
building during adaptive control of plants with 
unstructured model uncertainty. 
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The logic shown in figure (V-4) is the logic that was imple- 
mented for the model building example in the next chapter. It could 
be expanded to apply to wider classes of problems. The time needed 
to perform the model building precludes its use for unstable plants 
or where the speed of adaptation is critical. In the next subsec- 
tion, the ability to avoid problems resulting from on-line adapta- 
tion is discussed. 

Although not studied here, it may be possible to use the LMS 
algorithm in parallel with the IMLMS algorithm; the MS filter could 
be used to refine the estimated parameters of the modes that are 
being added by the IMLMS filter. Depending upon the application, it 
may be possible to restart the identification process after one com- 
plete model has been built using the refined estimates of the previ- 
ous model as initial estimates. This is possible because of the 
parallel structure of the plant, compensator, and IMLMS algorithm. 

E-3 Discussion of Additional Adaptive Problem Solving Logic 

It is desirable, especially for the case of model building, to 
have an. autonomous adaptive system that is capable of handling most 
adverse contingencies that may be encountered. A prime example is 
that of a space structure under construction. An adaptive control- 
ler would be required to satisfy certain pointing or damping 
requirements. Since it is envisioned that the entire control mis- 
sion would be satisfied by on-board computers, it will be necessary 
to expand the capabilities mentioned in the previous two subsections 
to include recovery from predictable (or even unpredictable) errors 
or divergences that may occur. 

In this subsection some potentially important factors in the 
adaptive design are considered. Although none of the problems 
described in this subsection were encountered during the simulations 
of the next chapter, some were encountered during other applications 
of this equivalent system identification scheme for adaptive con- 
trol. Other points are natural extensions that could conceivably be 
added. When a control system is fine tuned for use, it requires 
many iterations by a knowledgeable controls expert. Hence, it is 
envisioned that an ’’expert system" with knowledge to apply heuristic 
rules would be required to solve the problem. 
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This becomes especially feasible if multiple processors are 
used so that controller speed is not impacted while decisions are 
made. Looking at figures (1-8) and (1-9), a logical separation of 
functions would be to put the controller and recursive identifica- 
tion in one processor and the state space realization, compensator 
design, and adaptive logic rules in the other processor. 

An obvious protective mechanism would be the capability to turn 
off the compensator. Sometimes the overall scheme breaks down, 
either due to an inherent limitation or implementation error. This 
probably can initially be detected by a control input saturation or 
too large of an amplitude of output, the first indication of a 
divergence. Since this scheme, at this point, is intended for 
stable plants, the wisest thing would most likely be to stop con- 
trolling and restart the adaptive identification logic. A reason- 
able expection would be that the identification has failed (inaccu- 
rate) or the compensator model order is too small. 

A problem that immediately comes out of the analysis in section 
I-D is the scaling of the state variables and the control inputs 
with respect to each other and the step size factor, p. Otherwise, 
parameter divergence is likely. The codes for LMS and IMLMS algo- 
rithms of appendix II include normalization factors for the measure- 
ments and inputs. It may be necessary for the computer to adjust 
these factors during on-line, autonomous operations. This can 
normally be detected by 0 lms growing too large or by the para- 
meters of the LMS or IMLMS filter becoming unrealistically large. 

As it turns out, is not a very accurate measure of 

convergence. As indicated in section II-D, the error of prediction 
of the LMS or IMLMS algorithm is easily satisfied even prior to 
parameter convergence. So actually a LMS is an indication of 
divergence of the LMS or IMLMS algorithm. A better measure of con- 
vergence of the algorithm would be to check and see the rate-of- 
change of the parameters over a sufficiently long time. This is 
especially true for the control influence terms of the multivariable 
LMS or IMLMS algorithms, if a batch least squares method is not 
used for the numerator terms. 
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Another check for convergence of the LMS or IMLMS algorithms 
would be to evaluate the pole and zero locations of the discrete 
transfer function prior to realizing in state space form. It should 
be possible to isolate certain numerical difficulties prior to 
implementing a compensator. For example, sometimes a pole on the 
negative real axis in the z-domain may be identified during para- 
meter convergence problems. This corresponds to an imaginary fre- 
quency and usually results in an inaccurate model for the compensa- 
tor design. When the numerator terms are not converging, experience 
has shown that the discrete transfer function zeros frequently lie 
significantly outside the unit circle in the z-domain, corresponding 
to zeros deep in the right half plane of the s-plane. 

In the case of large space structures, it is expected that a 
significant period of time will be required to perform the model 
identification needed to construct a model for the control system 
design. In addition to using the control system actuators to excite 
the structure and making measurements from the control system sen- 
sors, it may be beneficial to utilize special purpose actuators and 
sensors to aid in the model identification process during this 
period of time. Application of disturbances other than through the 
actuators may be necessary to insure that the system has had all its 
modes excited. Auxiliary sensors could be used to distinguish 
between uncontrollable modes and sensor noise. The use of on-board 
heuristic logic would make the application of such practices 
feasible. 

Another task that could be turned over to an intelligent compu- 
ter system that is required to function autonomously, is the selec- 
tion of the designer's step size factor, p and the actual control 
(dither) inputs. These are the prime variables available to the 
controls engineer in this scheme for influencing overall system 
performance. Ideally, it would be advantageous to develop a learn- 
ing system that could vary p, u arnp and to enhance perfor- 

mance. Most likely a trial and error approach with systematic 
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variations would be required. So again, heuristic logic would need 
to be programmed based upon the designer's experience. 

Other factors could be considered. The fundamental issue is 
that a truly adaptive system will need some implementation of heur- 
istic rules to help prevent instabilities and reduce the risk of 
divergence. However, the proper balance between the vise of artifi- 
cial intelligence and control theory for solving real time control 
problems will need considerable research. 

Initially, fairly simple decisions need to be evaluated under 
the formalism of artificial intelligence as a means of verifying the 
approach. When the required decision process becomes highly com- 
plex, the approach can then be extended with confidence to enable an 
efficient search for the proper solution. 

In this research an attempt was made to include some heuristic 
rules to aid in the solving of off-nominal problems. What was 
actually implemented was described in the previous two subsections. 
In the next chapter some examples are presented where the overall 
scheme for adaptive control are applied to some generic problems. 
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Chapter VI 


ADAPTIVE CONTROL EXAMPLES 

A. RESTRUCTURABLE CONTROL EXAMPLE 

In this section a generic example of an unstable plant will be 
considered. The actuator will be assumed to break and an implemen- 
tation of the overall scheme proposed in this research will be used 
to restructure the control to achieve acceptable performance as 
described in section V-E.l. Notice that the terms that will be 
tracked by the adaptive filter are the numerator terms. It was 
shown in sections II-D and III-A that the adaptation of the control 
Influence terms are the slowest and most inaccurate of this 
approach. As a means of improving the adaptive controller perfor- 
mance, the batch least squares scheme of section II-E will also be 
used in conjunction with the recursive LMS filter. 

The linear, time-invariant plant which is representative of a 
highly unstable attitude hold mode for a helicopter in hover is 
given by 


( 6 . 1 ) 
( 6 . 2 ) 

where the system matrices are given in modal coordinates as 


with output 


x = Ax + Bu 


y ~ Hx 


'2 

A = 

.1 

-2" 

0 . 

(6.3) 

B = [2 

0] T 

(6.4) 

H = [0 

1] 

(6.5) 
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The open loop eigenvalues are s=l±j. 

The plant is clearly unstable and stabilization is required to 
maintain accurate regulation. A digital compensator is used, as 
described in chapter V. The estimator is 

x(k+l) = $x(lc) + Tu(k) - K[y(k)-Hx(k)] , (6.6) 


and the control law is 


u(k) = -Cx(k) 

+ Fy (k) , 

J com * 

(6.7) 

where, for a sample rate of 20 Hz 

» 


'1.10 

$ = 

..053 

-.105* 
.997 J 

(6.8) 

o 

• 

II 

t-H 

.00258 ] T . 

(6.9) 

The open loop poles in 

the z~domain are 

z=l .05±.05j 

(unstable). The constant compensator gain matrices. 

that were 

selected using the methods of chapter V, are 


K = [4.22 

1.20 ] T 

(6.10) 

C = [5.26 

15.4] 

(6.11) 

F = [ 

16.4] 

(6.12) 


The closed loop plant performance has eigenvalues z=.754±.16j. 
The apparent closed loop natural frequency of the plant is 6.68 
rad/ sec with a damping ratio of .48. 
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A simulation was performed using the approaches outlined in 
appendix I where the output is commanded by 


Y 

com 



sin(fit)>.767 
.707>sin(flt)>-.707 
sin(fitK -.707 


(6.13) 


with Y amp = .25 and ft=.25 rad/sec. 

At time t=25 seconds, half of a double "actuator” is assumed to 
break, resulting in only half the original control authority. That 
is, B changes to 


B = [l 0] T (6.14) 

The failure is sensed and the adaptive algorithms are given 50 
seconds to find a new representation of the plant using the 
recursive LMS algorithm. At time t=75 seconds a new controller is 
brought on-line in an effort to improve the control system 
performance. 

Time histories of the output and the commanded output are shown 
in figure (VI-1). The closed loop damping ratio goes from .78 to 
approximately .05 at 25 seconds. Stability is barely maintained and 
there is a large error in the commanded output. Although the system 
has not diverged, its performance is poor. 

After a preselected interval of 50 seconds has expired, the new 
compensator is brought on-line using the coefficients identified by 
the LMS algorithm. The system is stabilized with an apparent damp- 
ing ratio of approximately .45. However, there is still an error — a 
bias in the nonzero set point of the output command due to poor 
knowledge of B, a necessity in the computation of the feedforward 
gain matrix, F. 

If the LMS algorithm is given more time to find the numerator 
terms, it does a better job. In fact, a good strategy would be to 
update the compensator again after a certain interval in time. 
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OUTPUT TIME HISTORY 

NUMERATOR TERMS USING RECURSIVE LMS 



Figure VI-1 — Time history of output and commanded output for 

restructurable control example using recursive LMS 
algorithm. 


CONTROL INFLUENCE COEFFICIENT 



TIME, sec 

Figure VI-2 — Time history of actual control and estimated control 
influence coefficients using the recursive LMS and 
batch least square algorithms. 
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Figure (IV-2) shows time histories of the actual and estimated 
control influence term, . The estimate lags (as observed in 

section II-D); improvements are made only within a few samples after 
each nonzero pulse. 

Figure (VI-2) also shows a locus of estimated b]j from a 
single batch least squares estimate of the numerator terms as 
described in section II-E. It uses a minimum of the last 400 data 
points (20 seconds) or a maximum of the last 1000 data points (50 
seconds) . As can be seen, the batch least squares estimate for the 
numerator terms converges faster. At t=75 seconds the standard LMS 
has identified b^ with a 25 percent error, while the batch 
least squares approach has only a 5 percent error. 

A time history of the simulation using the batch least squares 
approach for identifying the numerator terms is shown as fig- 
ure (VI-3). When the new compensator is brought on-line, it 
stabilizes the system to nearly the same damping ratio as before the 
change in actuator characteristics and has a small commanded output 
error. 

The performance index for the controller is the RMS commanded 
output error (c.f. from eq. (5.35)). A plot of this parameter, 
normalized by y CO m> is shown as a time history in figure (VI-4). 
When the new compensator is implemented at t=75 seconds, the error 
is significantly reduced. However, if a goal of having o c under 
10 percent was set into the computer logic (see fig. (V-2)), then 
the adaptive scheme attempts to improve by implementing another 
compensator. It continues to do so until the goal is reached. The 
algorithm with just the LMS filter achieves this goal after t=120 
seconds. 

There are some weaknesses in the scheme. If b^ had become 
.25, say, then the closed loop system would have gone unstable. The 
LMS algorithm tends to diverge for unstable systems and the system 
may break before a new controller is brought on-line. However, with 
limits on the control input, u, the system may not be stabilizable 
anyway. Also, it take 50 seconds to bring an acceptable compensator 
on-line, which may not be fast enough for many applications. 
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OUTPUT TIME HISTORY 

LMS FILTER WITH BATCH LEAST 
SQUARES FOR NUMERATOR TERMS 



Figure VI-3 — Time history of output and commanded output for 

restructurable control example using recursive LMS 
algorithm in conjunction with batch least squares 
identification of the control influence terms. 


OUTPUT FOLLOWING PERFORMANCE OF RESTRUCTURABLE 

CONTROL EXAMPLE 



Figure VI-4 — Time history of normalized error of commanded output 
for restructurable control example. 
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B. MODEL BUILDING EXAMPLE 


In this example a low order model of a flexible spacecraft is 
considered. The rigid body mode is, of course, known and the flex- 
ible body modes are identified one at a time by the IMLMS algorithm 
and added to the control model. Although the plant is stable, 
acceptable control performance is possible only when an accurate 
model is used for the compensator design. 

Reference [67] shows that a sixth order model of the OSO-8 
spacecraft [68,69] is possible when near pole zero cancellations are 
taken into account. In fact, figure (VI-5) is a convenient repre- 
sentation of the model, which is similar to the study of flexibility 
upon control in reference [21]. 

Using the state space form of equations (6.1) and (6.2), the 
matrices for the system in figure (VI-5) are 


A = 


0 

m i 

0 

\ 

m„ 


1 

0 

0 

0 

0 

0 


0 

h 

m. 


■(k 1 +k 2 ) 

m 2 

0 

m„ 


0 

0 

1 

0 

0 

0 


0 

0 

0 

h 

m 2 

0 

-kj 

m„ 


0 

0 

0 

0 

1 

0 


(6.15) 


B = [0 0 

0 1 0 

of 

(6.16) 

SC 

II 

o’ 

o 

1 0 0 

0] 

(6.17) 

The 0S0-8 spacecraft 

had 

the following 

equivalent 


characteristics for the system of figure (VI-5) 


1 -1 

— = 13.96 sec id = 16.5 sec 

IUj 1 


id ^ = 40.64 sec * (6.18) 
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SIXTH ORDER SYSTEM EQUIVALENT TO REDUCED 
ORDER SPACECRAFT CONTROL PROBLEM 



Figure VI-5 — Diagram showing equivalent system description of 
of reduced order spacecraft control model. 
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For the convenience of the simulation, let time be measured in units 
of and mass in units of • The model parameters for the 
simulation are chosen to be 


m 


1 * *2 ' °i * 1 


(6.19) 



13.96 sec 


( 6 . 20 ) 


k x = .716 = 2.817 (6.21) 

Viscous damping corresponding to a damping ratio of .005 was 
added to the two modes. The simulation is performed using the 

techniques of appendix I at a sample rate of 20 Hz. The output is 
commanded through a pulsed wave given by (6.13) with y a mp = * 5 and 
&=.25 rad/ sec. 

A time history showing the output when a compensator was 
designed neglecting any knowledge of the flexible body modes is 
shown in figure (VI-6). The normalized RMS error of commanded 

output, o c , from equation (5.35) approached a steady state value 
of 18 percent. Additionally, the regulation about zero is 
apparently poor as a bias exists. 

The time history in figure (VI-7) shows the output if a compen- 
sator is designed including the first flexible mode. The RMS of 
commanded output error, o c , approaches a steady state value of 
.126 and the bias about zero is no longer a problem. While the 
closed loop system is effectively regulated, stable and damped, the 
remaining flexible mode is lightly damped and prominent in the out- 
put. Clearly, the response is unsatisfactory for any precision 
control requirement. 

Figure (VI-8) is a time history of the output when all three 
modes are included in the model to design the compensator. In this 
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OUTPUT VERSUS NORMALIZED TIME 

COMPENSATOR MODEL WITH 1 MODE OF EXACT DATA 



Figure VI-6 — Time history of output with rigid body mode modeled 
exactly for the suboptiraal compensator design. 


OUTPUT VERSUS NORMALIZED TIME 

COMPENSATOR MODEL WITH 2 MODES OF EXACT DATA 



Figure VI-7 — Time history of output with rigid body and first 
flexible body mode modeled exactly for the 
suboptimal compensator design. 
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case, a c approaches a steady state value of .083 which is repre- 
sentative of the good performance. There is no steady-state error 
and the modes are fairly well damped. However, use of the optimal 
control design strategies described in chapter V leads to notch 
filter designs for the compensators [70,71], Other methods of 
designing the control system may in fact lead to superior perfor- 

mance, especially given inaccuracies in knowledge about the 
dynamics. 

Designing control systems for spacecraft can be difficult [73] 
because ground-based testing generally leads to inaccurate represen- 
tations of spacecraft performance. Furthermore, large space struc- 
tures may require months to construct. Although certain pointing 
and damping tasks may be required during this deployment stage, it 
would be difficult to model the configuration at each possible, 
intermediate step. For these reasons, it is desirable to have 
learning or adaptive control schemes similar to the one suggested by 
this research. 

The IMLMS algorithm was used to build a model for this 

example. The rigid body dynamics are known and the output was com- 
manded as in figure (VI-6) for 50 units of normalized time. At that 
time, the coefficients from the IMLMS algorithm are used to add the 
first flexible mode to the model used to design the compensator. 
The next 50 time units are shown in figure (VI-9). The normalized 
a c has been reduced from 18.2 percent to 13.1 percent, which is 
similar to the 12.6 percent of the compensator with exact knowledge 
for two modes only. Figures (VI-7) and (VI-9) look fairly similar. 
The steady-state error has been removed, but the undamped oscilla- 
tion of the third mode still exists. 

After 100 time units, the IMLMS algorithm is used to add the 

second flexible mode to the model. A time history of the ensuing 
motion is shown as figure (VI-10). The performance is nearly as 

good as the exact case shown in figure (VI-8). The RMS error of 
commanded output, o c , is equal to 10.1 percent versus 8.3 percent 
for the exact case. The output is well regulated and the modes are 
fairly well damped. 
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OUTPUT VERSUS NORMALIZED TIME 

COMPENSATOR MODEL WITH 3 MODES OF EXACT DATA 



Figure VI-8 — Time history of output with plant modeled exactly 
for suboptimal compensator design. 


OUTPUT VERSUS NORMALIZED TIME 

COMPENSATOR MODELWITH 2 MODES OF IMLMS DATA 



Figure VI-9 — Time history of output with rigid body mode modeled 
exactly and first flexible mode obtained from IMLMS 
filter for suboptimal compensator design. 
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OUTPUT VERSUS NORMALIZED TIME 

COMPENSATOR MODEL WITH 3 MODES OF IMLMS DATA 



Figure VI-10 — Time history of output with rigid body mode modeled 
exactly and both flexible modes obtained from IMLMS 
filter for suboptimal compensator design. 


121 



The compensator that results from the optimal control design is 
of the same order as its system model and contains a notch filter. 
This filter is highly sensitive to identification errors. Other 
design methods (e.g. [70,71]) may not be as sensitive to the para- 
meter estimation accuracy of the flexible body modes, but may give 
poorer performance. 

Simulations were also made using the LMS algorithm where all 
modes are identified at the same time. An improvement over no know- 
ledge of the flexible body modes was possible using the LMS algo- 
rithm to identify and add these modes to the compensator model. 
However, small parameter errors translate into large dynamics errors 
as the order of the system to be identified increases. In fact, 
when the order of the system was assumed to be four, the LMS filter 
identified a negative real z-plane pole which is quite inaccurate. 

The results of these simulations are summarized in table VI-1 - 
The IMLMS algorithm used in conjunction with the model building 
scheme improves the accuracy and damping of the commanded response. 
One reason such good estimation accuracy is possible is because some 
prior knowledge about the flexible body modes is used in selecting 
the step-size scale factors of the IMLMS algorithm. The denominator 
of the discrete transfer function for a single damped mode can be 
written as [21] 


d(z) = z^ - 2e a ^(cosbT)z + e (6.22) 

T is the time between samples, a is the real part of the 

s-domain roots and b is the imaginary part of the s-domain roots. 
Since space structures are lightly damped, the real part of any pair 
of unaugmented roots will be approximately zero. So, if the coef- 
ficients of the LMS or IMLMS filters correspond to the following 

denominator equation, 


d(z) 


2 

z 


V 


- a r 


(6.23) 


the parameters can be equated as 
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Table 6.1 — Performance of adapted and nonadapted compensators in controlling the 
output of a reduced order model of the 0S0-8 spacecraft. 


SOURCE OF MODEL 

EXACT 

I MI 

MS 

LM 

S 

NO. MODES IN MODEL 

1 

2 

3 

2 

3 

2 

3 

EIGENVALUES OF 
MODE 
1 

0,0 

0,0 

0,0 

0,0 

0,0 

NA 2 

0,0 

EIGENVALUES OF 
MODE 
2 

NA 1 

-.005 ± j 

-.005 ± j 

-.006 ± 1.04j 

-.006 ± 1 .04 j 

-.071 ± 1 .19 j 

.006 ± 1.23J 

EIGENVALUES OF 
MODE 
3 

NA 1 

NA 1 

-.012 ± 2 .52 j 

NA 1 

-.013 ± 2.52j 

NA 1 

-.19 ± 2.3j 

NORMALIZED ERROR 
OF COMMAND 
FOLLOWING 

.182 

.126 

.083 

.131 

.101 

.169 

.157 


NA 1 Not available, mode is not included in truncated compensator model 

NA 2 Not available, root found by LMS has imaginary frequency (negative real root on z-plane) 




cij = 2(cosbT) 


(6.24) 


« 2 - -1 (6.25) 

Since 0.2 is approximately -1 with good accuracy (especially for 
high sample rates), it can be preset to -1 and the appropriate step 
size factor nulled. There is no need to identify a term that is 
known a priori. This, in effect, reduces the number of degrees of 
freedom, thereby improving the speed of adaption and the accuracy of 
convergence. 

Even better accuracy for the model building adaptive controller 
using the IMLMS could be achieved by allowing more time for conver- 
gence and varying the input signals (especially pulse frequency) in 
a systematic way. A logical expectation is that minutes or possibly 
hours could be taken by a learning system to identify a model of 
high order to perform a control mission in space on a structure 
under modification or construction. 

A verification of the MIMO capabilities of the algorithm was 
performed by adding a colocated actuator-sensor pair at mass m 3 . 
The number of free parameters in the denominator polynomial coeffi- 
cients remained unchanged so the identification accuracy remained 
approximately the same. Adequate pulsing of both controls for both 
outputs required more time to identify the coefficients of the 
numerator polynomial. Equivalent accuracy for the numerator terms 
was obtained in twice the time as used for the SISO system. As pre- 
viously illustrated, the full model built by using the IMLMS algo- 
rithm resulted in fairly similar output as using the exact model in 
the compensator design. The inclusion of the additional sensor and 
actuator improved the damping and yielded a c T s of .038 and .051 
for the exact model and IMLMS sequentially identified models, 
respectively. Hence, no problems were encountered extending the 
scheme to MIMO systems. 
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ADAPTIVE CONTROL IN THE PRESENCE OF SENSOR NOISE 

Many of the adaptive control algorithms that have been studied 
diverge in the presence of observation noise (c.f. chapter I). A 
small bias (constant or oscillatory) in the sensor output produces 
divergence. This instability is avoided by separating the system 
identification and the controls tasks. Only occasionally is the 
information allowed to flow to the control block. Hence, the adap- 
tive control scheme of this research will not exhibit any worse per- 
formance than that from using a time-invariant, suboptimal con- 
troller. In this example, the adaptive control scheme proposed in 
this research is used to augment sensor performance. 

Assume for this example a linear, time-invariant plant which is 
known exactly and has the following familiar form 

x_ = Ax + Bu (6.26) 

with output 

y = Hx + s^ + s g sin(s (0 t) (6.27) 


where s^, is the sensor bias, s a is the sensor noise amplitude 
and s u is the sensor noise frequency. The bias and noise are, of 
course, uncontrollable by u. The system matrices were chosen from 
an unstable Dutch roll mode of an airplane and are given by: 



(6.28) 


B - [0 2] T (6.29) 

H = [2 0) (6.30) 


The open loop eigenvalues are s=-.l±j, damping ratio of .1 and a 


natural frequency of 1 rad/sec. 
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A digital compensator of the form of equations (6.6) and (6.7) 
was designed to control the plant output. The equivalent discrete 
system, for a sample rate of 20 Hz, has system matrices 



' .999 

.0491 


$ = 

.-.049 

.989 J 

(6.31) 

r = [.002 

.099] T 

(6.32) 


The gain matrices for the digital compensator, designed using 
the methods of chapter V, are 

K = [.519 .361 ] T (6.33) 

C = [-55.8 -7 .40] (6.34) 


F = [28.2] 


(6.35) 


The closed loop eigenvalues are z=.586 + .289 j . The output is 
commanded to follow equation (6*13) with yamp = *^^ anc ^ rad/ 

sec. 

In this example, an extreme case is chosen with s^ equal to 
.05, s a equal to .05, and s^ equal to 3. Figure (VI-11) shows 
the actual and commanded output for s a =01 • The corresponding 
sensor output is plotted in figure (VI-12). The control effort 
expended to reduce the oscillation in the sensor output results in 
exciting an oscillation in the plant output, which is clearly not 
desirable. 

After 100 seconds of recursive identification using the IMLMS 
algorithm, an extra mode is added to the compensator model. The 
sensor noise frequency was identified within 2 percent and the 
control influence terms were very small, less than 10"^, making 
the mode uncontrollable for the tolerance levels selected in the 
software design package [51]. The result of applying the new model 
to the system is shown in figures (VI-13) and (VI-14). 
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PLANT OUTPUT 

WITH UNMODELED SENSOR NOISE 



Figure VI-11 — Time history of plant output with oscillatory 
sensor noise. 


SENSOR OUTPUT WITH UNMODELED SENSOR NOISE 



Figure VI-12 — Time history of sensor output with oscillatory 
sensor noise. 
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PLANT OUTPUT 

WITH MODELED UNCONTROLLABLE MODE 



TIME, t sec 


Figure VI-13 — Time history of plant output with oscillatory 
sensor noise. A mode has been added to the 
compenstator model from the IMLMS filter. 


SENSOR OUTPUT WITH MODELED UNCONTROLLABLE MODE 



TIME, t sec 


Figure VI-14 — Time history of sensor output with oscillatory 
sensor noise. A mode has been added to the 
compensator model from the IMLMS filter. 
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The plant output in figure (VI-13) remains virtually unchanged, 
so the performance is not impacted by knowledge of the sensor 
noise. However, the sensor output is significantly better as shown 
in figure (VI— 14). This is a result of using output in the cost 
function for finding the control law; this causes the extra mode to 
be included as if it were a physical mode of the plant. While the 
perceived performance is enhanced, the actual performance is poor, 
indicating the need to identify and minimize sensor errors. 

The problem can be solved if the oscillatory mode can be 
classified as sensor noise. This is accomplished by testing the 
identified mode to see if it satisfies some threshhold of control- 
lability. Equation (5.11) is then modified such that no ouput from 
the uncontrollable mode appears in the y used in the cost function 
for computing the controller feedback gains. This makes y equiva- 
lent to the plant output instead of the sensor output. The results 
from this experiment are shown in figures (VI-15) and (VI-16). The 
plant output has relatively good performance and the sensor output 
has reverted to oscillations. In effect, the uncontrollable but 
detectable mode is subtracted from the total sensor output to esti- 
mate the plant mode. This made it possible to observe and control 
the actual plant output in spite of the observation noise. However, 
the problem would still exist for cases where prior knowledge does 
not allow one to distinguish between sensor noise and an uncontrol- 
lable mode of the plant. 

The influence of sensor bias is shown in figures (VI-17) and 
(VI-18). The plant output is well damped, stable and follows the 
commanded output well, except that it is offset from the desired 
output by sjj. The sensor output is also offset slightly, but to a 
smaller extent. The bias was not included in the control model and 
the threshholds that were chosen for the adaptive logic did not 
result in the addition of another mode when used in parallel with 
the IMLMS algorithm. If good following performance is desired, the 
sensor bias needs to be small. 

A main observation from this example is that the present adap- 
tive control scheme does not diverge in the presence of uncontrol- 
lable modes. Thus, if some heuristic decision v logic can be 



SENSOR OUTPUT WITH SENSOR NOISE MODELED 



Figure VI-15 — Time history of plant output when uncontrollable 
mode is modeled as sensor noise. 


PLANT OUTPUT WITH MODELED SENSOR NOISE 



Figure VI-16 — Time history of sensor output when uncontrollable 
mode is modeled as sensor noise. 
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PLANT OUTPUT WITH SENSOR BIAS 
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Figure VI-17 — Time history of plant output with an unmodeled 
sensor bias. 


SENSOR OUTPUT WITH SENSOR BIAS 



Figure VI-18 — Time history of sensor output with an unmodeled 
sensor bias. 
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implemented to distinguish between an uncontrollable mode of the 
plant and observation noise, it is possible to use this adaptation 
scheme to improve plant performance in the presence of such sensor 
noise. 


132 



CHAPTER VII 


CONCLUDING REMARKS AND RECOMMENDATIONS 
A. SUMMARY 

A scheme for the adaptive control of multi-input/multi-output 
(MIMO) plants is proposed that utilizes heuristic logic. Control 
law restructuring takes place intermittently after satisfying a 
number of convergence tests instead of continuously. The plant is 
constructed by starting with a simple model and identifying one 
additional mode at a time until the added mode makes no substantial 
improvement in closed-loop performance. The recursive identifica- 
tion scheme which allows the adaptation of part of the assumed model 
while holding the remaining part fixed is a modification of the LMS 
(Least Mean Square) algorithm. It is termed the incremental mode 
LMS ( IMLMS) algorithm. 

On-line identification is performed using an extension of the 
LMS algorithm of Widrow and Hoff to MIMO systems. This is a recur- 
sive identifier having the same form as the extended Kalman filter 
(EKF), but it uses constant instead of time-varying gains. Hence, 
it does not require the propagation of a large covariance matrix, 
and it avoids the "oblivious filter" problem that occurs with the 
EKF when the variance of the estimated parameters becomes small. 
Excitation for identification is performed using rectangular pulse 
out put- command s . 

Only the coefficients of the denominator polynomial of the MIMO 
z-transfer function are identified using the IMLMS algorithm. The 
coefficients of the numerator polynomials are identified using a 
batch least squares algorithm over the last several output com- 
mands. Control law synthesis is done by LQG methods which require 
the formation of a MIMO state space realization from the z-transfer 
function identified by the IMLMS algorithm. An efficient regression 
algorithm is presented to do this and is illustrated by a second 
order example with two inputs and two outputs. 
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Three adaptive control examples were presented: (1) a control 
restructuring example where a partial failure of the actuator 
resulted in a change in actuator gain; (2) a model building example 
of a spacecraft with a rigid body mode and two flexible body vibra- 
tion modes; and, (3) a sensor noise example where uncontrollable 
oscillatory noise appears in the sensor signal, but not in the plant 
output being controlled. 

As in previous adaptive control schemes, there is no way of 
discerning whether detectable but uncontrollable modes appear in the 
actual plant output or only in the sensor output. If this determin- 
ation can be made by some other means, then an appropriate model can 
be formed. 

The proposed scheme for adaptive control is not suitable for 
situations where rapid adaptation is required; in particular, it 
will not handle plants with significant instabilities. 

B. CONCLUDING REMARKS 

Several continuously adapting control algorithms have been 
proven to be globally stable. It has been shown, however, that the 
assumptions required for proving stability are overly restrictive. 
In fact, divergence is likely when operated with sensor noise, with 
biases in controls or measurements, or with unmodelled dynamics 
[17-20] present. The goal of this research was to extend a simple 
recursive estimation algorithm to a form which will minimize the 
impact of these real world considerations. This is accomplished 
with the proposed adaptive control scheme, but the added complexity 
prevents proving global stability. 

The LMS algorithm was selected as the baseline parameter esti- 
mator because of its low computational burden and good convergence 
characteristics in the presence of process noise, sensor noise and 
unmodelled dynamics. The algorithm was enhanced by the systematic 
inclusion of MIMO systems. A prime contribution of this research is 
the development of the IMLMS algorithm which provides a systematic 
means for including unmodelled dynamics. The IMLMS algorithm is a 
recursive identification scheme that estimates the parameters of 
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part of the assumed model while holding the remaining part of the 
model constant. This provided the basis for an on-line structure 
identification scheme using heuristic rules to build models. 

Although linear equivalent systems are sought as models, plants 
with small or locally insignificant nonlinearities will be readily 
controlled by the proposed adaptive control scheme. Plants with 
large nonlinearities which prevent the ready representation as a 
linear system will not be easily controlled by this or any adaptive 
controller now being considered in the research literature. Excep- 
tions would be cases where complex plants are modelled with good 
physical understanding and have only a few time varying or uncertain 
parameters. In this research, it is assumed that very little prior 
knowledge about the plant is available, necessitating model learning 
and building in real time. The approach in this research is useful 
for plants which can be effectively linearized, but have a wide 
range of uncertain modal characteristics or are slowly changing. 

The model used in the compensator for controlling the outputs 
was updated only after satisfying a number of convergence tests and 
passing performance threshholds. Ry preventing the controller gains 
from being updated at each sample, the divergence problems that 
plague continuously adaptive control systems are avoided. Unfortu- 
nately, this approach of waiting to evaluate performance and conver- 
gence properties on-line prevents direct analysis. Without the 
ability to directly analyze the overall scheme, proofs for global 
stability are impossible to obtain. However, this research does 
provide practical solutions to the generic types of problems that 
adaptive control algorithms suffer from. Although mathematical 
rigor is replaced with engineering judgement, this scheme of using 
heuristic logic points to a class of learning controllers which may 
have high reliability for a wide range of applications. 

A final concern of the adaptive control scheme proposed in this 
research is the added complexity of utilizing heuristic logic to 
avoid divergence. Software reliability is already a difficult 
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problem and an area of intense research. The use of simple hypothe- 
sis testing may make the issue of program quality assurance a sig- 
nificant hindrance in its use for other than highly experimental or 
back-up applications of adaptive control. Additionally, the proper 
balance between control theory and machine intelligence concepts 
needs to be considered in terms of problem solving capability, com- 
puter requirements, problem application and software reliability. 

C. RECOMMENDATIONS FOR FURTHER STUDY 

The prime recommendations for developing enhanced adaptive con- 
trol methodologies are those involved with making adaptive control a 
practical alternative to robust control. As a first step, the 
results reported in this document should be verified through hard- 
ware implementation. A dual microprocessor control experiment for a 
lightly-damped, flexible beam using the algorithms of this report is 
being considered for the experimental apparatus described in refer- 
ence [73]. It will provide a benchmark for comparing with other 
real-time adaptive control experiments. Trade-offs between recur- 
sive identification and batch processing will be possible, plus the 
computational efficiency of performing the compensator design can be 
evaluated. It could also be an effective check of the approach of 
using heuristic adaptive logic for dealing with practical problems 
as they arise. It may be more cost effective to perform hybrid 
(digital and analog) simulations than trying to perform accurate, 
high-order, digital simulations. 

One advantage of hardware experimentation is that the LMS adap- 
tive filter was originally developed for analog implementation. It 
may be possible to find a similar way to implement the IMLMS filter 
in a hybrid control system. Analog implementation of the recursive 
identification techniques could provide faster adaptation, more 
accurate parameter estimation and less sensitivity to the numerical 
convergence, especially at high sample rates. However, the use of 
analog adaption will introduce sampling errors into the parameter 
estimates during digitizing. 
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Lattice filters [74] are a potentially attractive alternative 
to the LMS class of adaptive algorithms considered in this 
research. It has been shown that an estimate for the system order 
can be recursively included in the identification scheme [75] , 
making rapid model identification possible. If n is the model 
order, the lattice filter algorithm has approximately lOn more com- 
putations per cycle than the LMS algorithm [74]. It does, however, 
provide the capability to iterate the model order based upon hypo- 
thesis testing at frequent, recursive intervals. Research shows 
that it can be useful for the identification of space structures 
[76] , but its use for cases with substantial damping has not been 
verified. Furthermore, direct use of the lattice filter output for 
compensator design is still of concern and needs further research. 
Another shortcoming is that the lattice filter cannot explicitly 
include any prior knowledge of the plant. In spite of these limita- 
tions, the use of lattice filters with their model order iteration 
capability and low computational burden does look to be promising 
for future adaptive control research. 

The use of quadratic cost techniques may not be ideal for 
on-line compensator design. An alternative design technique should 
be considered for the observer and controller implementation. It 

may be possible to utilize computationally efficient and robust 
algorithms that take advantage of the modal forms preserved from the 
recursive identification schemes. In addition, the feedforward con- 
troller design algorithm should be extended to handle a different 
number of control inputs from output measurements [64] and to be 
able to follow other than step inputs with arbitrary placement of 
the zeros [21]. 

An area of active research that all adaptive control approaches 
would benefit from is the continued development of a method for 
determining optimal control inputs for system identification [43- 
45]. Implementation of a full dual control scheme — where control 
power expenditures for improving identification, as well as for 



performance requirements, are explicitly included in the cost func- 
tion that is optimized — appears to be too computationally burden- 
some at this time. However, an inexact technique may be possible 
whereby suboptimal inputs are computed or suggested by the adaptive 
strategies implemented in the computer by designers. 

In this vein of "smart" controllers, whole new avenues need to 
be explored. It may be possible to couple an "expert" system from 
the science of artificial intelligence utilizing rules and knowledge 
from experienced controls engineers to help build an on-line adap- 
tive controller. An adaptive control scheme is envisioned which 
could be successfully applied to a wide variety of problems while 
requiring a minimum of problem dependent fine-tuning and still main- 
tain a high probability of success. However, a logical next step 
would be merging the science and rules under development for artifi- 
cial intelligence with the formal analysis frequently applied to 
control system design. This may eventually result in a flexible and 
fast "universal controller" which could be used to improve control 
system performance for many types of problems. 
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APPENDIX I 


DESCRIPTION OF SIMULATION EXPERIMENTS 

A, ADVANCED CONTINUOUS SIMULATION LANGUAGE (ACSL) 

The computer simulations were performed using the CDC Cyber 
network at Langley Research Center. ACSL (Advanced Continuous Simu- 
lation Language) [77] was the prime tool used in these studies. It 
provides a facility for performing integration of nonlinear equa- 
tions of motion and for obtaining a wide variety of outputs from 
either interactive or batch use. ACSL interprets the set of com- 
mands given it, which look like FORTRAN, and writes a FORTRAN pro- 
gram for the user. Hence, it is relatively easy to couple FORTRAN 
subroutines to the ACSL program. 

The simulations shown in this thesis were performed using a 
second order Runge-Kutta integration scheme with a step size at 
least as small as the sample time of the discrete compensator and 
identification schemes. The program was segmented so the required 
field length at execution would be small enough to allow interactive 
execution. Graphic output was obtained from Tektronix copiers 
during interactive simulation and Varian plotters for batch process- 
ing. 

Subroutines were added to the ACSL programs to provide discrete 
sampling with zero order hold, the compensator and control law 
implementation, performance evaluation and the heuristic adaptive 
decison making. A brief description of these routines is given in 
table 1-1. Program flexibility was maintained making it relatively 
easy to interchange models and adjust parameters during the debug- 
ging stages of the adaptive control law development. Input decks 
minimized the need to recompile the entire simulation during problem 
formulation and evaluation. 


Table 1-1 


Subroutines used to augment ACSL for 
simulating adaptive control schemes. 


SUBROUTINE 

DESCRIPTION 

ADMODE 

Adds modes to compensator model 
during model building process. 

COMMAG 

Selects input pulses for 
identification. 

COMND 

Distributes pulses to appropriate 
control or output command. 

CONTROL 

Computes control input from feedback 
and feedforward law. 

ERREVAL 

Accumulates figure-of-merit 
computation for adaptive control. 

FORM 

Assembles dynamics matrix from modal 
information. 

IMLMS 

Performs IMLMS algorithm computation, 
(see app. II) 

KFILT 

Computes state estimates using 
asymptotic KBF gains. 

LMS 

Performs LMS algorithm computation 
with MIMO extensions (see app. II) 

SAMPCK 

Code for implementing discrete 
sampling logic. 

SETUP 

Routine for initializing plant and 
compensator model selection. 
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B. OPTIMUM REGULATOR ALGORITHMS FOR THE CONTROL OF 
LINEAR SYSTEMS (ORACLS) 


The code used for performing the matrix operations used by ACSL 
and the compensator designs as described by chapter V is ORACLS 
[51] • ORACLS is a library of subroutines for synthesizing control 
logic for multivariable state space systems. ORACLS also provides a 
set of routines for matrix operations and manipulations. 

ORACLS was coupled to ACSL making on-line compensator design 
possible. However, a number of routines had to be added to custom- 
ize the library to the present problem. In addition, several analy- 
sis routines were added to aid in debugging the adaptive control- 
ler. Descriptions of the subroutines that were used to augment the 
ORACLS library are given in table 1-2. 

Arrays for ORACLS are stored as a single subscripted array 
stacked by columns in combination with an integer pointer array. 
For example, if a 3 X 2 matrix A is to be used by a routine from 
ORACLS, it is saved as array A(l) of length 6, where the first three 
elements are the first column of A and the last three elements are 
the second column. In addition, an integer array NA(J) of length 2 
is also used. The first value is the number of rows and the second 
value is the number of columns. 
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Table 1-2 — Subroutines used to augment ORACLS to facili- 
tate on-line implementation for adaptive 
control schemes. 


SUBROUTINE 

DESCRIPTION 

COMPEN 

Computes gain matrices for asymptotic, 
suboptimal compensators. 

CONVRT 

Transforms roots between s-plane and 
z-plane . 

CTRLAW 

Computes control law for asymptotic, 
discrete, suboptimal regulator. 

ELEMENT, 
GRAB, ETC. 

Several matrix manipulation routines 
to faciliate the use of ORACLS. 

INIT 

Initializes common blocks and local 
variables for ORACLS and ACSL. 

KFILTD 

Computes gain matrices for asymptotic, 
discrete Kalman-Bucy filter. 

MATINV 

Finds general matrix inverse including 
pseudo inverses for non-square matrix. 

PHIGEN 

Computes discrete equivalents for 
linear system matrices. 

POLZER 

Finds poles and zeros of linear system 
on both s-plane and z-plane. 

SETPT 

Computes the feedforward gain matrix 
for non-zero set point of regulator. 
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APPENDIX II 


COMPUTER CODES FOR OUTPUT ERROR IDENTIFICATION 

This appendix contains FORTRAN listings of the LMS and IMLMS 
recursive, adaptive filters. The subroutines are called at each 
discrete time step to update the predicted output and the parameter 
estimates. YP is the predicted output. The storage arrays for com- 
puting the ARMA predictions are XS and US and are updated each time 
step. XNORM and UNORM are used to normalize the inputs to the 
filter as a means of scaling the data to near unity as suggested in 
chapters II and III. The parameters being updated are ALPHA and 
BETA, which are the coefficients of denominator and numerator poly- 
nomials, respectively. 
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Subroutine listing of FORTRAN code for LMS Adaptive Filter. 


SUBROUTINE LMS ( IORD , XMU , ALPHA , BETA , Y , NY , YNORM , U , NU 


2 ,UNORM,XS,US,YP,ERR,IOUT) 

a************************************************************************ 

* SUBROUTINE LMS IS THE LMS ADAPTIVE FILTER USED AT EACH * 

* TIME STEP. IT IS USED FOR ADAPTIVE ESTIMATION * 

* AND PARAMETER IDENTIFICATION. * 

************************************************************************* 

« 

* 

* INPUTS: 


* 

* 

« 

* 

* 

* 

* 

* 

* 

* 

* 

* 


IORD IS THE ORDER OF THE DENOMINATOR 
XMU IS THE STEP SIZE FACTOR 
ALPHA ARE THE DENOMINATOR COEFFICIENTS 
BETA ARE THE NUMERATOR COEFFICIENTS 
Y IS THE MEASUREMENT OF THE OUTPUT 
U IS THE INPUT (CONTROL) 

XS IS A STORAGE ARRAY OF LAST Y’S 
US IS A STORAGE ARRAY OF LAST U'S 
XNORM AND UNORM ARE NORMALIZING FACTORS 
NY AND NU ARE THE DIMENSIONS OF Y AND U 
IOUT=1 FOR OUTPUT 


* 

# 

* 

* 

* 

* 

» 

» 

* 

* 

* 

* 


* OUTPUTS: # 

* 

* ALPHA & BETA ARE UPDATED EACH TIME STEP * 

* XS AND US ARE UPDATED EACH TIME STEP * 

* YP IS THE PREDICTED OUTPUT FOR THE CURRENT TIME * 

* ERR IS THE ERROR OF PREDICTION FOR THE CURRENT TIME * 

a************************************************************************ 


REAL ALPHA(I) ,BETA(2,2,2) ,Y( 1 ) ,U( 1 ) ,XS(2 ,2) ,US(2,2) 
REAL YP( 1 ) ,ERR( 1 ) 

INTEGER NY ( 1 ) ,NU(1) 

IMEAS=NY( 1 ) 

ICTRL=NU ( 1 ) 


C 

C COMPUTE PREDICTION 
C 

DO 200 L=1 .IMEAS 
XSTO=0 . 

USTO=0 . 

DO 190 1=1, IORD 
XSTO=XSTO+ALPHA(I)*XS(I ,L) 

190 CONTINUE 

DO 195 1 = 1, IORD 
DO 195 J=1 .ICTRL 
USTO=USTO+BETA (I,J,L)*US(I,J) 


195 CONTINUE 

YP(L)=USTO+XSTO 
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200 CONTINUE 


COMPUTE ERR & NEW WTS 


* 


DO 250 L=1,IMEAS 
ERR(L)=Y (L)-YP(L) 

250 CONTINUE 

DO 300 1=1 , IORD 
DO 300 L=1 ,IMEAS 

ALPHA( I)=ALPHA(I)+2 .*XMU*XS(I ,L)*ERR(L) 

300 CONTINUE 

DO 305 J=1 ,ICTRL 
DO 305 1=1, IORD 
DO 305 L=1 ,IMEAS 

BETA(I,J,L)=BETA(I,J,L)+2.*XMU*US(I,J)*ERR(L) 

305 CONTINUE 
C 

C UPDATE MEAS & CTRL STATES 
C 

DO 320 L=1 ,IMEAS 
DO 310 11=2, IORD, 1 
I=I0RD+2-II 
XS(I ,L)=XS(I-1 ,L) 

310 CONTINUE 

XS(1 ,L)=Y(L)/YNORM 
320 CONTINUE 

DO 330 J=1 ,ICTRL 
DO 325 11=2, IORD, 1 
I=I0RD+2-II 
US(I,J)=US(I-1 ,J) 

325 CONTINUE 

US(1 ,J)=U(J)/UNORM 
330 CONTINUE 


C 

C 

C 


DEBUGGING OUTPUT 


IF(IOUT.EQ.O) GO TO 500 
WRITEC6 ,410) IORD , IMEAS , ICTRL 
410 FORMAT (/5X" LMS ALGORITHM OUTPUT 

1 ,I5/10X,"N0. OF MEAS . = " ,I5/10X,"NO. 
WRITE (6 ,415) (ALPHA(I) ,1=1 .IORD) 

4 1 5 F0RMAKT20 ," ALPHA : "5 (T20 , 5F 1 2 . 5/) ) 
WRITE(6 ,420) 

420 FORMAT ( 1 0X , " BETA : " ) 

DO 430 J=1 , ICTRL 


"/10X,"SYSTEM ORDER 
OF CTRLS.= ” ,15) 


DO 430 L=1, IMEAS 

WRITE(6 ,425) J ,L, (BETA(I , J ,L) ,1=1 ,IORD) 

425 FORMAT( 10X,"CTRL= " ,15 ,5X ,"IMEAS= ",I5,5(T50,5F12.5.)> 


430 CONTINUE 

WRITE(6 ,440) 

440 FORMAT ( / 1 0X , "STATES : ” ) 

DO 450 L=1, IMEAS 

WRITE(6 ,445) L,(XS(I,L) ,1=1 ,IORD) 
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445 FORMAT( 10X,"MEAS= " , 15 ,5(T30 .5F12.5, ) ) 

450 CONTINUE 

WRITE(6 ,455) 

455 FORMAT ( / 1 OX , " CONTROL STATES:") 

DO 460 J=1 .ICTRL 

WRITE(6 ,458) J,(US(I,J) ,1=1 ,IORD) 

458 FORMATC 10X, "CONTROL= ",I5 ,5(T30 ,5F12 .5 ,) ) 
460 CONTINUE 
500 RETURN 
END 
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Subroutine listing of FORTRAN code for IMLMS Adaptive Filter 


SUBROUTINE IMLMSC IDEN , INUM , XMUD ,XMUN , ALPHAF , BETA , IDAD , ALPHAI 
2 ,Y ,NY ,YNORM ,U ,NU ,UNORM,XS ,US ,YP ,ERR , IOUT) 

******* ************************* ****** ************ ****************** 


* SUBROUTINE IMLMS IS THE INCREMENTAL MODE LMS * 

* ALGORITHM APPLIED AS A RECURSIVE IDENTIFIER * 

* ON PART OF THE MODEL. * 

******************************************************************** 

* * 

* INPUTS: * 

* * 


* IDEN IS ORDER OF FIXED PART OF DENOMINATOR * 

* INUM IS ORDER OF NUMERATOR * 

* XMUD IS ARRAY OF DENOMINATOR STEP SIZE FACTORS * 

* XMUN IS ARRAY OF NUMERATOR STEP SIZE FACTORS * 

* ALPHAF ARE COEFFICENTS OF FIXED PART OF DENOM. * 

* BETA ARE NUMERATOR COEFFICIENTS * 

* IDAD IS ORDER OF UNKNOWN PART OF DENOM. * 

* ALPHAI ARE THE COEFFICENTS OF INCREMENTAL PART * 


* Y IS OUTPUT MEASUREMENT * 

* US IS CONTROL INPUT * 

* NU AND NY ARE DIMENSIONS OF U & N * 

* UNORM AND YNORM ARE NORMALIZING FACTORS * 

* XS AND US ARE STORAGE ARRAYS * 

* IOUT=1 FOR OUTPUT * 

* * 

* OUTPUTS: » 

* * 

* ALPHAI & BETA ARE UPDATED EACH TIME STEP * 

* XS & US ARE UPDATED EACH TIME STEP • 

* YP IS THE PREDICTED OUTPUT AT CURRENT TIME * 

* ERR IS THE ERROR OF PREDICTION * 


******************************************************************** 
REAL ALPHAF(I) ,BETA(2,2,2) ,Y( 1 ) ,U( 1 ) , XS(2 ,2) ,US(2 ,2) 

REAL YP ( 1 ) , ERR( 1 ) .ALPHAI ( 1 ) ,XMUD( 1 ) ,XMUN( 1 ) 

INTEGER NY Cl) ,NU( 1 ) 

IMEAS=NY( 1 ) 

ICTRL=NU(1) 

NEWD=IDEN+IDAD 

IF(IOUT.EQ.I) WRITE(6 ,410) IDEN , IMEAS , ICTRL 
IDELAY = 1 

COMPUTE PREDICTION 

DO 200 L=1, IMEAS 
XST01=0. 

XST02=0. 

XST03=0. 

UST01 =0 . 
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non no 


DO 190 1=1 , IDEN 

XST01 =XST01 +ALPHAFC I ) *XS ( I ,L) 

DO 190 K=1 , IDAD 

XST03 =XST03-ALPH AF ( I ) * ALPHAI (K ) «XS ( I+K , L ) 

190 CONTINUE 

DO 192 K=1,IDAD 
XST02=XST02+ALPHAI (K ) *XS(K ,L) 

192 CONTINUE 

DO 195 1 = 1 , INUM 

DO 195 J = 1 .ICTRL 

USTOI =UST01+BETA( I ,J,L)*US(I,J) 

195 CONTINUE 

YP ( L ) =UST01 +XSTO 1 +XST02+XST03 
200 CONTINUE 

COMPUTE ERR & NEW WTS 

DO 250 L=1 ,IMEAS 
ERR(L)=Y(L)-YP(L) 

250 CONTINUE 

COMPUTE NEW PARAMETERS 

IF(XS(NEWD,1).EQ.O) GOTO 309 
DO 300 K=1 , IDAD 
DO 300 L=1 ,IMEAS 
T1 =0 . 

F=2.*XMUD(K)*ERR(L) 

DO 299 1=1 1 IDEN 
T1 =T1-ALPHAF( I)*XS(I+K ,L) 

299 CONTINUE 
UPDATE=F* (XS(K ,L)+T1 ) 

ALPHAI (K) = ALPHAI (K)+UPDATE 
IF(IOUT.EQ.O) GOTO 300 

WRITEC6 ,40) K.ALPHAI(K) ,XS(K,L) ,T1 .UPDATE 
40 FORMAT( 10X"K, ALPHAI, XS, T1 ,UPDATE== " ,15 , 4F15 .6) 

300 CONTINUE 
IF(IOUT.EQ.O) GOTO 301 

WRITE ( 6 , 42 ) YP ( 1 ) , Y ( 1 ) , ERR ( 1 ) ,XST01 , XST02 , XST03 , USTO 1 
42 FORMAT ( 1 OX"YP , Y ,ERR="3F12 .6 ,5X, "XSTO( 1 ,2,3)=".3F12.6, 
1 5X,"USTO=",F12.6) 

C 

c 

301 DO 305 J=1 .ICTRL 
DO 305 1 = 1 , INUM 
DO 305 L=1 ,IMEAS 
F=2.*XMUN(I)*ERR(L) 

BETA( I , J ,L)=BETA(I , J ,L)+F*US(I , J) 

305 CONTINUE 
C 

C UPDATE MEAS & CTRL STATES 
C 

309 DO 320 L=1,IMEAS 
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DO 310 11=2, NEWD, 1 
I=NEWD+2-II 
XS(I,L)=XS(I-1 ,L) 

310 CONTINUE 

XS(1 ,L)=YP(L)/YNORM 
320 CONTINUE 

DO 330 J=1 ,ICTRL 
IF(INUM.LT.2) GOTO 326 
DO 325 11=2 , INUM , 1 
I=INUM+2-II 
US(I,J)=US(I-1 ,J) 

325 CONTINUE 

326 US(1 ,J)=U(J)/UNORM 
330 CONTINUE 

DEBUGGING OUTPUT 

IF(IOUT.EQ.O) GO TO 500 

410 F0RMAT(/5X" LMS ALGORITHM OUTPUT " /I OX ."SYSTEM ORDER 
1 , 15/1 OX ,"N0. OF MEAS.= " ,I5/10X, M NO. OF CTRLS.= ",I5) 

WRITE (6 ,415) ( ALPHAI(I) ,1=1 , IDAD) 

4 1 5 FORMATC 1 0X , " ALPHAI : "5 (T20 ,5F12 .5/) ) 

WRITEC6 ,420) 

420 FORMATC 1 0X, "BETA:") 

DO 430 J=1 .ICTRL 
DO 430 L=1 .IMEAS 

WRITE (6 ,425) J ,L , ( BETA ( I , J ,L) , 1=1 ,INUM) 

425 FORMATC 1 0X, "CTRL= " ,15 ,5X,"IMEAS= " ,15 ,5(T50 ,5F12 .5 , ) ) 

430 CONTINUE 

WRITEC6 ,440) 

440 FORMATC 1 0X, "STATES:") 

DO 450 L=1 .IMEAS 

WRITEC6 ,445) L,(XS(I,L) ,1=1 ,NEWD) 

445 FORMATC 1 0X, "MEAS= " ,15 ,5(T30 ,5F12 .5, ) ) 

450 CONTINUE 

WRITEC6 ,455) 

455 FORMATC 1 0X, "CONTROL STATES:") 

DO 460 J = 1 .ICTRL 

WRITEC6 ,458) J,(US(I,J) ,1=1 ,INUM) 

458 FORMATC 10X,"C0NTR0L= " ,15 ,5 (T30 ,5F12 .5 , ) ) 

460 CONTINUE 
500 RETURN 
END 
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